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FOREWORD 


This report presents results of the expansion and improvement of the 
FORMA system for response and load analysis. The acronym FORMA stands 
for FORTRAN Matrix Analysis. The study, performed from 16 May 1975 
through 17 May 1976 was conducted by the Analytical Mechanics Department, 
Martin Marietta Corporation, Denver Division, under the contract NAS8- 
31376. The program was administered by the National Aeronautics and 
Space Administration, George C. Marshall Space Flight Center, Huntsville, 
Alabama under the direction of Dr. John R. Admire, Structural Dynamics 
Division, Systems Dynamics Laboratory. 

This report is published in seven volumes; 

Volume I - Programming Manual, 

Volume IIA - Listings, Dense FORMA Subroutines, 

Volume IIB - Listings, Sparse FORMA Subroutines, 

Volume IIC - Listings, Finite Element FORMA Subroutines, 

Volume IIIA - Explanations, Dense FORMA Subroutines, 

Volume IIIB - Explanations, Sparse FORMA Subroutines, and 
Volume IIIC - Explanations, Finite Element FORMA Subroutines. 
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ABSTRACT 


This report presents techniques for the solution of structural 
dynamic systems on an electronic digital computer using FORMA (FOR TRAN 
Matrix Analysis). 

FORMA is a library of subroutines coded in FORTRAN IV for the effi- 
cient solution of structural dynamics problems. These subroutines are 
in the form of building blocks that can be put together to solve a large 
variety of structural dynamics problems. The obvious advantage of the 
building block approach is that programming and checkout time are limi- 
ted to that required for putting the blocks together in the proper order 

The FORMA method has advantageous features such as: 

1. subroutines in the library have been used extensively for many 
years and as a result are well checked out and debugged; 

2. method will work on any computer with a FORTRAN IV compiler; 

3. incorporation of new subroutines is no problem; 

4. basic FORTRAN statements may be used to give extreme flexi- 
bility in writing a program. 

Two programming techniques are used in FORMA: dense and sparse. 
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I. INTRODUCTION 

A listing of the source deck of each sparse FORMA subroutine is 
given in this volume to remove the "black-box” aura of the subroutines 
so that the analyst may better understand the detailed operations of 
each subroutine. 

The format of a sparse matrix on a utility tape is giv^n in Chapter 
II. 

The FORTRAN IV programming language is used in all sparse FORMA 
subroutines. 
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II. FORMAT OF SPARSE (Y ) MATRIX ON UTILITY TAPS ''DISK) 


/ 


end of record 


Matrix 

Partition 1 j 

Partition 1 

Partition 2 

Partition 2 

etc j 

Header 

Header } 

Numbers 

Header 

Numbers 

1 


— 

(2 records) 


(2 records) 
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MATRIX HEADER : 

HROW 

NCOL 

NPART 

NNZA 

IFORD 

KV 

ISHAPE 



Number of rows in matrix 
dumber of columns in matrix 

Number of partitions of matrix on tape (disk) 

Number of non-ssrces in matrix 
Indicator for ordered matrix 

Dimension size of work vector when matrix was formed 
Shape indicator (DIAG. LOWER, UPPER, WHOLE) 

extra 


PARTITION 1 HEADER: 


NNZP 

LFELP 

LLEL? 

0 

0 

0 

0 

0 

0 

0 


Number of nonzeross in 
Location (Row, Column) 
Location (Row, Column) 


partition 

of first element in partition 
of last element in partition 


extra 


PARTIT ION 1 NUMBERS : 

(LV(I) , I ~ 1, KNZP) 

Location (Row, Column) 

(V(I), 1=1, NNZP) 

t 

Element Value 
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III. SUBROUTINE LISTINGS 

The subroutines are given in alphabetical order with numbers coming 
before letters. 



oooooooooooooonooooo 


YAA 


SUBROUTINE YAA ( ALPHA ,NU1 A ,NUTZ ,V, LV,KV ,NUT1 ,NUT 2 ) 

DIMENSION VU).LV(l) 

DATA NITfNOT/5 *6/ 

SCALAR ALPHA TIMES SPARSE MATRIX A. (ALPHA * A - Z). 

CALLS FORMA SUBROUTINES YIN ,YINI ,YLORD ,YNOZER,YOUT ,YOUTI , 

YPART »ZZBOMB • 

DEVELOPED PY R A PHILIPPUS. JANUARY 1970. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

ALPHA = SCALAR THAT MULTIPLIES MATRIX A. 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

NUT1 - LOGICAL NUMBER OF UTILITY TAPE. 

NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

NERROR EXPLANATION 
1 = SIZE LIMITATION EXCEEDED. 


REWIND NUTA 
REWIND NUTZ 

CALL YINI (NUTAyLVt 1*10) 
CALL YCUTI (NUTZ, LV, 1*10) 
NPAPT: LV(3) 

DO 10 J=1,NPART 

CALL YINI (NUT A »LV ,1,10) 

CALL YOUTI (NUT Z ,LV ,1,10) 

NNZ=LV(1) 


C 


C 


C 


IF (NNZ.GT.KV) GO TO 999 
IF (NNZ.GT.O) GO TO 3 
CALL YOUTI (NUTZ ,LV ,1,1) 


CALL YOUT 
GO TO 10 
3 CALL YINI 
CALL YIN 


(NUTZ, 


V , 


1 , 1 ) 


(NUT A ,LV , 1 ,NNZ ) 
(NUT A ,V, 1 ,NNZ ) 


DO 5 1=1, NNZ 
5 V(1)=ALPHA*V(I ) 


CALL YOUTI (NUTZ, LV,1, NNZ) 
CALL YOUT (NUTZ ,V, 1 »NNZ ) 
10 CONTINUF 


CALL YNOZER (NUTZ ,V,LV ,KV,NUT1 ) 

CALL YLORD (NUTZ ,V,LV,KV,NUT1 ,NUT2) 
RETURN 


NERR0R=1 


999 CALL ZZBOMB (3HYA A 
END 


,NERROR ) 
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SUBROUTINE YAABB (ALP HA » NUTA ,B ETA , NUTB , NUT 2 »V «IV« KV»NUT1 ,NUT2 ) 
SUBROUTINE YAAPB (ALP HA , NUTA , BETA, NUTB , NUTZ ,V , LV ,KV,NUT1 ,NUT2 ) 
SUBROUTINE YAABB (ALP HA »NUTA ,B ETA, NUTB, NUTZ ,V * LV ,KV,NUT1 ,NUT2 ) 
DIMENSION V(1),LV(1), MHE AD ( 1C ) 

DATA NIT, NOT/5, 6/ 

C 

C MATRIX SUMMATION FOR SPARSE MATRICES. ALPHA * A + BETA * B = Z. 

C CALLS FORMA SUBROUTINES XLCRD ,YIN ,Y.’NI » YLORD ,YNOZER,YCUT 

C YOUTI » YPART ,YSYMLH, YSYMUH t YZERO. 

C DEVELOPED BY R A PHILIPPUS. MAY 1969. 

C LAST REVISION BY RL WOHLEN FOR NASA. MAY 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C ALPHA = SCALAR THAT MULTIPLIES MATRIX A. 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

C BETA = SCALAR THAT MULTIPLIES MATRIX P. 

C NUTB = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX B IS STORED. 

C NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 * LOGICAL NUMBER OF UTILITY TAPE. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C GET (A) HEADER INFORMATION. 

REWIND NUTA 

I REWIND NUTB 

I CALL YINI (NUTA »MHEAD, 1 »1C ) 

NR A = MHE AD ( 1 ) 

NCA = MHE AD (2) 

NNZ A = MHEAD (4 ) 

IASHAP = MHEAD( 7 I 

C GET (B) HEADER INFORMATION. 

CALL YINI (NUTB, MHEADv 1,10) 

NRB = MHFAD(l) 

NCP = MHEAD(2) 

NNZP = MHE AD (A ) 

IPSHAP =• MHE AD ( 7 ) 

C 

C ALLOW FOR DIFFERENT SIZE (A) AND (B). 

NRZ=NRA 

NCZ-NCA 


IF 

(NRB . 

GT. NR A) NR Z 

=NPR 





IF 

(NCB . 

GT. NCA) NCZ 

-wee 





IF 

(NNZ A 

.GT. 0 .AND. 

NNZB 

.GT. 0) 

GO 

TO 

4 

IF 

(NNZ A 

.EO. 0 .AND. 

NNZB 

.GT. C) 

GO 

TO 

15 

IF 

(NNZ A 

• G T • 0 • AN D • 

NNZB 

.EQ. C) 

GO 

TO 

15 


CALL YZERO (NUT Z ,NRZ ,NCZ ) 

RETURN 

C MAKE (A) AND (P) SAMF SHAPE. 

4 IF (IASHAP.Nf .IBSHAP) GO TO 5 
IZSHAP^IASHAP 

GO TO 15 

5 IZSHAP=5HWHOLE 

IF (lASHAP.cQ. 5H WHOLE .OR. I ASHAP.E0.4HDI AG) GO TO 1C 
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IF (IASHAP. E0.5HLCWER ) CALL YSYMUH (NUTA, V ,LV .KV.NUT1 ,NUT2 ) 
IF ( IASHAP. EQ.5HUPPER ) CALL YSYMLH CVUTA,V,LV,KV,NUT1,NUT2) 
10 IF (IBSHAP. EC. 5HWHOLE .OR. IBSHAP. E0.4HDIAG) GO TO 15 

IF (IBSHAP. EQ.EHlfWER ) CALL YSYMUH (NUTB , V ,LV ,KV ,NUT1 ,NUT2 ) 
IF ( IBSHAP. E0.5HUPPER ) CALL YSYMLH (NUTB * V»LV »KV, NUTI »NUT2 ) 
C MAKE CERTAIN ELEMENTS ARE ORDERED • 

15 CALL YLORD (NUTA,V,LV,KV,NUT1 ,NUT2 ) 

CALL YLORD (NUTB ,V, L V,KV, NUTI ,NUT2 ) 

C 

LZS-1 
1=0 
1A=0 
IB=0 
NREC=0 
REWIND NUTA 
REWIND NUT6 
REWIND NUTZ 
C 

CALL Y1NI (NUTA *MHE AD *1,10) 

NPARTA r MHEAD ( 3 ) 

NNZA = MHEAD (A ) 

IASHAP = MHEAD (7) 

CALL YINI (NUTB .MHEAD, 1,10) 

NPARTB = MHEAD (3 ) 

NNZB = MHEAD (A) 

IBSHAP = MHEAD (7) 

MHEAD(l) = NRZ 
MHEAD (2 ) = NCZ 

IF (NNZA.GT.O .AND. NNZB.GT.O) GO TO 35 
C 

IF (NNZA.GT.O) GO TO 25 

MHEAD ( 3 ) = NPARTB 

MHEAD(A-) = NNZB 

MHEAD(7) = IBSHAP 

CALL YOUTI (NUTZ, MHEAD, 1,10) 

C 

DO 20 J*l, NPARTB 

CALL YINI (NUTB ,MHE AD, 1,10) 

CALL YOUTI (NUTZ, MHEAD, 1,10) 

NMULT = MHEAD(l) 

CALL YINI (NUTB ,LV, 1 ,MHEAD( 1 ) ) 

CALL YIN (NUTB ,V * 1 , MHEAD ( 1 ) ) 

DO 18 IMULT=1, NMULT 
18 V(IMULT) = BETA*V (IMULT) 

CALL YOUTI (NUTZ, IV, 1,MHFAD(1)) 

20 CALL YOUT (NUTZ, V,1,MHEAD(1)) 

RETURN 

C 

25 MHEAD (3 ) = NPARTA 
MHEAD(A) * NNZA 
MHEAD(7) = IASHAP 
CALL YOUTI (NUTZ .MHEAD, 1,10) 

DO 30 J=l, NPARTA 

CALL YINI (NUT A ,MHE ADt 1 , 10 ) 

CALL YOUTI (NUTZ ,MHE AD, 1 , 10 ) 
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NMULT = MHEADd ) 

CALL YINI (NUTAtLV,ltMHEAD(l)) 

CALL YIN CNUTA ,Vtl ,MHEADU) ) 

00 28 IMULT=1*NMULT 
28 V1TMULY) = ALPHA*V(IMULT ) 

CALL YOIITI (NIJTZ »LV* 1 ,MHEAD 11)) 

30 CALL YOUT (NUTZ »V, 1 *MHEAr ( 1 ) ) 

RETURN 

35 NRFAD=N» > ART A+NPARTB 
NZMAX=KNZ A+NNZE 
MHEAD (3) = NRFAD 
MHEAD(A) = NZMAX 
MHEA0C6) = C 
MHEAD (7 ) = 1ZSHAP 
CALL YOUTI (NUTZ »MHE AD* 1 » 10 ) 

CALL YINI (NUTA,MHEAD*1»10) 

NNZPA = MHEADd) 

LFELPA = *«hE AD (21 
LLELPA = MHE ADO I 
CALL YINI (NUT8 »MHE AD » 1 » 10 ) 

NNZPB = MHEADU) 

LFELPE = MHE AD (2 ) 

LLELPE - MHEAD ( 3 ) 

READ A PARTITION OF A AND MULTIPLY IT BY ALPHA 
40 1=1+1 
IA=IA+1 

LZE = LZS-1 +NN7 0 A 

CALL YINI (NUTA,LV,LZS*LZE ) 

CALL YIN (NUTA,V»LZS*LZE) 

IF (IA.OF .NPaRIA) GO TO 42 
CALL YINI (NUTA,MHEAD,1 ,10) 

NNZPA = MHf. (I) 

LFELPA = MHEAD ( 2 ) 

LLFLPA = MHE ADO) 

42 IF (ALPHA. FT.. I. I GO TO 50 
DO 45 J=LZS,LZE 
45 V(J)=ALPHA*V(J) 

50 IF (I.GT.l) GO TO 65 
LZS=LZE+1 

READ A PARTITION OF P AND MULTIPLY IT PY BETA 
55 1=1+1 
IB=IP+1 

LZE=LZ5-1 +NNZPB 

CALL YINI (NOTE *LVfLZS »LZF ) 

CALL YIN (NllTB *V* LZS * LZE ) 

IF (IB.GE.NPARTP ) GO TO 57 
CALL YINI (NUTB»MHEAD,1,10) 

MNZPE = MHFAD(1 ) 

LFFLPB = MHFAD ( 2 ) 

LLELPB = MHE AD ( 3 ) 

57 IF (BETA.E0.1.) GO TO 65 
DO 60 J=LZS »LZ t 
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60 V( J)=BETA*V( J) 

65 LAF=LZ 5-1 
IZ=LZS 
C 

DO 85 IP=1,LAE 

70 IF CLV(IP |-LV(IZ I J 85,80,75 
75 IZ=IZ*1 

IF (IZ.GT.LZE) GO TO 90 
GO TO 70 

80 V(IP)=V(IP)*V(IZ> 

V(IZ)=C. 

IZ*IZ+I 

IF (IZ.CT.LZE) CO TO 90 
85 CONTINUE 
C 

90 NNZW=0 

c 

DO 95 IZ=I »LZF 
IF (V(IZ).FQ.O.) GO TO 95 
NNZW=NNZW ♦! 

V (NNZW ) =V ( IZ ) 

LVINNZW )=LV I IZ ) 

95 CONTINUE 
C 

CALL XLCRL (V,LV,1,NNZW) 

LZE=NNZW 

WRITl A PARTITION OF Z 
100 MAX=KV/4 

IF (LZE.LT.KV/4*3 .A*'D. I.LT.NREAD) GO TO 130 
IF (MAX.GT.LZE) MAX = LZE 

105 IF (LV(MAX) . LT.LFELPB .OR. I B . EO.NP APTB > GO 110 
MAX=MAX-1 

GO TO 105 

110 IF (LV(MAX).LT .LFELPA .OR. I A.EQ.NPARTA ) GO TO 115 
MAX=MAX— 1 
GO TO 110 

115 IF (MAX.EO.D GO TO 120 
MHEAD( 1 ) = MAX 
MHFAO (?) = LV( 1 ) 

MHEAD(3) = LV(MAX) 

MHEAD(IO) = 0 

CALL YOUTI (NUTZ ,MHEAD,1 ,10 ) 

CALL YOUTI (NUT Z,LV,1,MAX) 

CALL YOUT (NUTZ ,V, 1 , MAX) 

NREC=NREC+1 

DETERMINE WHETHER AND WHAT TO READ OR WRITE 
120 IF (MAX .EQ.LZE .AND. I.EO.NREAO) GO TO 135 
K>MAX 

MOVE -LZE— MAX 
DO 125 J=1,MCVE 
K=K ♦? 

V ( J ) =V( K ) 

1?5 L V( J ) = L V( X ) 
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LZF=MPVE 

IF (T.EC.NREAD) CP TO 100 
130 LZS=LZF+1 
MIK=LFFLPA 

IF (MIN.GT.LFELPB .AND. IE .LT.NPARTB) MIN=LFE LPB 
IF (MIN.EO.LFELPA .AND. LZE+NNZPA.GT.KV) G" TO ICO 
IF (MIN.EC'.LFEIPB .AND. LZE+NNZPB.GT.KV) GO TO 100 
IF (MIN.EO.LFELPA .AND. I A.LT. NPARTA ) GO TO 40 

GO TO 55 

135 IF fNRFC.EO.NPEAD) GO TO 140 
DC 138 J=1,1C 
138 MHEAD(J) = 0 

CALL YOUTI (NUTZ, MHEAD, 1,10) 

CALL YOUTI (NUTZ ,MHE AD* 1 * 2 ) 

CALL YOUTI (NUTZ ,MHEAD,1 *2) 

NREC=NRFC + 1 
GC TO 135 

140 CALL YNOZER (NUTZ ,V,L V,KV,NUT1 1 
RETURN 
END 
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SUBROUTINE YASSEM (NUT A, IR Z, JC Z.NUTZ ,V, LV, KV,NUT 1 ,NUT2,NUT3 ) 
DIMENSION V(1),LV(1), MHE AD (101 
DATA NIT»N0T/5 *6/ 

C 

C SPARSF MATRIX ASSEMBLY. (MATRIX A INTO MATRIX Z). 

C BE SURE MATRIX 2 IS DEFINED BEFORE CALLING THIS SUBROUTINE. FOR 
C EXAMPLE, CALL YZERO TO CLEAR MATRIX Z. 

C CALLS FORMA SUBROUTINES YIN ,YINI , YLCRD *YOUT ,YCUTI ,YPART , 

C YSYMLH,YSYMUH,ZZBOME. 

C DEVELOPED BY P. A PHILIPPUS. JANUARY I97C. 

C LAST REVISION eY WA EENFIELD. MARCH 197fc. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 
C IRZ = ROW NUMBER IN MATRIX 2 OF FIRST ROW OF MATRIX A. 

C JCZ = COLUMN NUMBER IN MATRIX Z OF FIRST COLUMN OF MATRIX A. 

C NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

C V = VFCTOR WORK SPACE. 

C LV = VECTOP WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT2 = LOGICAL NUMBEP OF UTILITY T APE • 

C NUT3 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NEPROR EXPLANATION 

C 1 = MATRIX A EXCEEDS MATRIX Z. 

? ~ SIZE LIMITATION EXCEEDED. 

^ 3 = SIZE LIMITATION EXCEEDED. 

C 

REWIND NUTA 

CALL YINI (NUTA ,MHE AD ,1,10 

NR A = MHFADd ) 

NCA = MHE AD ( 2 ) 

IASMAP = MHEADC7 ) 

ILIMIT=NRA+IRZ-J 
JLIMIT=NC A+JCZ-1 
REWIND NUTZ 

CALL YINI (NUTZ , MHE AD, 1 ,10) 

NRZ = MHEAD(l) 

NCZ = MHE AO ( ? ) 

1ZSHAP = MHF AD ( 7 ) 

IZSAVF = MHFAD ( 7 ) 

NERR0R=1 

IF ( IL1MIT.GT.NRZ .OR. JLIMIT.GT.NCZ ) GO TO 999 
IF (IASHAP.EQ.5HL0WER ) CALL YSYMUH (NUTA,V ,LV ,KV,NUT1 ,NUT2) 

IF ( IZSHAP.EO.SHLOWFR ) CALL YSYMUH (NUTZ ,V ,LV ,KV , NUT 1 ,NUT2 ) 

IF (IASHAP.FC. 5H UPPER * CALL YSYMLH (NUTA,V ,LV,KV,NUT1 ,NUT2 ) 

IF (IZSHAP.E0.5HUPPER ) i ALL YSYMLH (NUTZ , V ,LV , KV, NUT1 ,NUT2 ) 
IZSHAP=5HWH0LE 

IF (IZSAVF. E0.4HDIAG .AND. I ASHAP . E0.4HDI AG .AND. IRZ. EQ. JCZ) 

* 1ZSHAP=4HDIAG 
REWIND NUI A 
REWIND NUT) 

CALL YINI (NUT A, MHE AD, 1, 1C) 

1ASHAP = MHE AD (7 ) 
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MHEAD(l) = NPZ 

MHEAD(2) = NCZ 

MHEAD(7 ) = IZSHAP 

CALL YOUTI CNUT1 , MHEAD, 1,10) 

NPARTA = MHEAD ( 3 ) 

NNZA = MHEAD (4) 

IF INNZA.GT.O) GC TO 2 
DO 1 1=1,10 
VCI) = 0. 

1 MHE ADC I ) = C 

CALL YOUTI (NUT1 ,MHE AD, 1 ,10) 
CALL YOUTI (NUT 1 ,MHEAO, 1,21 
CALL YOU 7 (NUT1, V,I,2) 
NPARTA=0 
GO TO 12 

2 LVADD=100CCC* ( IRZ— I ) ■*■ JCZ-1 
C 

DO 10 1=1, NPARTA 

CALL YINI (NUTA»LV* 1,10) 

IF (LVCD.GT.KV) GC TO 999 
LVC2 ) = L V( 2 )+LVADD 
LV(3)=LV(3)+LVADD 
CALL YOUTI CNUT1 ,LV, 1,10) 
NNZ=LV( I ) 

CALL YINI (NUT A *LV, 1 »NNZ ) 

, CALL YIN (NUTA,V,1 ,NNZ) 

DO «, J=I,NNZ 
5 LV(J)=LV( J)+LVADD 

CALL YOUTI (NUT1 ,LV, 1 »NNZ ) 

10 CALL YOUT (NUT1 ,V, 1 ,NN2 ) 

C 

12 REWIND NUTZ 
REWIND NUT? 

CALL YINI (NUTZ,MHEAD,1,IC) 
NPARTZ = MHE AD (3) 

NNZZ = MHEAD (4 ) 

MHEAD ( 7 ) = I2SHAP 

CALL YOUTI (NUT?,MHEAD,1 ,10) 

IF (NNZZ.LE.O) GO TO 30 


DO 25 1 = 1 ,NPA r TZ 
CALL YINI (NUTZ, LV, 1,10) 

IF (LV(l).GT.KV) GO TO 999 
CALL YOUTI (NUT2,LV, I ,10) 
NNZ=LV(1 ) 

CALL YINI (NUTZ ,LV, 1 ,NNZ ) 
CALL YIN (NUT Z ,V, 1 , NNZ ) 

DO 20 J=1 »NNZ 
IZ=LV( J )/lCCCOO 
IF (IZ.LT.IRZ) GO TO 2G 
IF (IZ.GT.ILIM1T ) CO TO 2G 
JZ=LV( J ) -ICC 000* I Z 


NERR0R=2 


NERR0R=3 
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IF OZ.LT.JCZ) CO TO 20 
IF UZ.GT.JLIMIT) GO TO 20 
V( J)=P. 

20 C0NT1IUE 

CAL’. YOUTI (NUT2,LV,1 ,NNZ) 

25 CALI YOUT (NUT2,V,1 ,NNZ) 

CELL YNOZER (NUT2 ,V,LV,KV»NUT3 ) 
GO TO 35 
3G REWIND NUTZ 
MHEA^Ill = NRZ 
MHEA'f2) = NCZ 
HHFA (31 = NPAPTA 
RHEA »U> = NNZ A 
MHFA.H5) = C 
MHE AO ( 6 ) = 0 
WH£AD(7) = IZSHAP 
CALL YOUTI (NUTZ »MHEADt 1*10) 

60 Tt’ 45 
35 REWIND NUT 2 

CALL YINI (L'UT2,MHEAD,1,10) 
NPART2 = MHEAD(3) 
mil - HHFAD(4) 

NP AF T= NPARTA-*NPAP TZ 

NNZ-NNZA+NNZZ 

REMIND NUTZ 

MH-A0I3) = NPART 

MHEADtA) = NNZ 

MHEAO ( 5 ) = C- 

MHEAD (*• 1 = 0 

MHE AD ( 7 ) = IZSHAP 

CALL YOUTI (NUTZ,MHEAD,1,10) 

DO 40 I =1 »NDAPTZ 

CALL YINI (NUT 2 »LV » 1 » 1 0 ) 

CALL YOUTI (NUTZ »LV* 1 *10) 
NNZ=LV( II 

CALL YINI (NUI 2 tLV* 1 *NN2 I 
CALL YIN (NUT2,V,1 ,NNZ) 

CALL YOUT . (NUTZ *LV« 1 *NNZ I 
40 CALL Yr-.iT (Nl'TZ ,V, 1 , NNZ ) 

45 REWMD NUT1 

CALL YINI (NUT1 *MHE AD * 1C »1 0 ) 

-...fcADnoi = o 

DC 50 1 = 1 .NPART A 

CALL YINI (NUT1 ,LV* 1 *10) 

CALL V’UiI (NUTZ »LV 1 1*10) 

NN? = LV( 1 ) 

CA_L YINI (NUT 1 t LVf 1 tNNZ ) 

CALL YIN (NUT1,V,1 ,NNZ) 

CALL YOUTI (NUTZ ,LV, 1 , NNZ ) 

50 CALL YOUT (NUTZ *V, I ,NNZ ) 
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CALL YLORD (NUTZ t V t L V»KV, NUT1 ,NUT2 > 
RETURN 

999 CALL ZZBOMB (6HYASSEM »NFPROR ) 

END 
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SUBROUTINE YBSL3A ( NU TU,NUTD *NUTB,NUTZ, V,LUV, KV.NUT1 ,NUT2 ) 
DIMENSION V(l), LUVC1), IH(IO) 

C 

C SPARSF BACK SOLUTION OF (U**I )*(-D-)*(U)*( Z) = « P I - 
C (U) IS A BANDED UPPER TRIANGULAR MATRIX WITH ONES ON THE DIAGONAL. 

C (-D-) IS A DIAGONAL MATRIX. IB)* AND THUS (Z) IS A MATRIX OF FULL 
C COLUMNS. ASSUMING (U)*(Z)=IY) AND (-D-)* (Y )= I G) GIVES ( U**T) * ( G) = ( B ) 
C WHICH ARE EASILY SOLVED FOR (G) t (Y), AND (Z). 

C (B)*(G)*(Y)«(Z) IN 1ST THIRD OF V. (U) IN 2ND THIRD OF V. 

C (-D-) IN 3RD THIRD OF V. 

C CALLS FORMA SUBROUTINES YIN ,YINI , YLORD *YOUT ,YOUTI ,YPART * 

C YTRANS. 

C DEVELOPED BY P L WOHLEN AND R A PHILIPPUS. MARCH 1972. 

C LAST REVISION BY R A PHILIPPUS. MARCH 1975. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTU - LOGICAL NUMBER OF UTILITY TAPE WITH MATRIX t*. 

C NUTD = LOGICAL NUMBER OF UTILITY TAPE WITH MATRIX D. 

C NUTB = LOGICAL NUMPER OF UTILITY TAPE WITH MATRIX B. 

C NUTZ = LOGICAL NUMBER OF UTILITY TAPE WITH CALCULATED MATRIX Z. 

C V = VECTOR WORK SPACE. 

C LUV = VECTOP WO p K SPACE. 

C KV = DIMENSION SIZE OF V*LUV IN CALLING PROGRAM. 

G NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT 2 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

E CONVERT B FROM SPARSE NOTATION TO FULL COLUMN NOTATION. 

• CALL YTRANS (NUTB ,NUT1 * V*LUV,KV,NUT2*NUTZ ) 

REWIND NUT1 
REWIND NUT2 

CALL YINI (NUT1 *IH* 1*10) 

NCB = I HC 1 ) 

NRB = IH( 2 ) 

NPART = I H ( 3 ) 

NNZP = IH(4) 

KVMN = KV— NR P 
KVMN02 = KVMN/2 
NCG = KVMN02/NPB 
NGB = (NCP-D/NCG+l 
DO 1 1=4,7 

1 IH ( I ) = 0 
IH(I ) = NPP 
IH ( 2 ) = NCB 
I H ( 3 ) = NGB 

CAlL YOUTI (NUT?, IN, 1 ,10) 

NNZRB = 1 

IE (NCB. LT. NCG) NCG=NCB 

LBS -- KV/A+1 

LBSM1 = LBS-1 

LBE = LBSM1 +NCG*NPB 

DO 2 I=LBS,LBE 

2 V(I) = C. 

JF = 1 

JL = NCG 

DO B 1=1 ,NP API 
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CALL YINI (NUT1 *IH* 1*10) 

NNZPB = IH(1) 

CALL YINI (NUT1 ,LUV , 1 ,NNZPB ) 

CALL YIN (NUT1, V,l, NNZPB) 

00 8 LB = 1 »NNZPP 

JB - lUV(LB)/1GOOOO 

IB = LUV( LB )-l0000C*JB 

IF (JB.LE.JL) GO TO 4 

IH ( ? = NCG 

IH ( 2 ) - 0 

IH(3) = C 

CALL VOUTI (NUT2,IH,1,10) 

CALL YOl/T (NUT2, V,LBS » LB E ) 

DO 3 J=LPS,LBE 

3 VCJ) = C. 

JF = JL+1 

JL = JF+NCG-1 
IF (JL.GT.NCB) JL=NCB 
NCG = JL-JF+1 
LBE = LBSM1+NCG*NRB 

4 L = ( JB-JF)*NRB+IB 
LBSM1L = LBSM1+L 
V(LESMIL) = V ( LB ) 

IF (I.LT.NPART .OR. LB. LT. NNZPB) GO T E 
IH(1) = NCG 
IH(2) = 0 
IH(3) = 0 

CALL YOUTI (NUT 2 ,IH,1,10) 

CALL YOUT ( NUT 2 * V,LES » LBE ) 

8 CONTINUE 
C 

C V(1 THRU (KV-NJ/2) CONTAINS B,G*Y,Z COLUMNS OF A GROUP. 

C V((KV-N)/2+l THRU KV-N) CONTAINS COLUMNS OF U (FROM DIAGONAL UP TO 
C TOP NON-ZERO) OF A GROUP. 

C V(KV-N+1 THRU KV) CONTAINS D. 

C LUV( 1 ) » 1=1 »N IS NUMBER OF ELEMENTS IN COLUMN I. 


REWIND NLTU 

CALL YINI (NUTU,IH, 1 ,10) 
N = I H ( 1 ) 


NGU = I H ( 3 ) 

LSU = (KV-N)/2 
LSD = KV-N+1 
CALL YINI 
REWIND N'UTD 
CALL YIN 
RFWIND NUT2 
CALL YINI 
NGB = IH( 3) 
REWIND NUT1 
CALL YOUTI 


+ 1 

(NUTU ,LUV , 1 ,N ) 

(NUTD »V»L SD , LSD+N-1 ) 
(NUT2,IH, 1,10) 


(NUT1 ,IH, 1,10 


DO 89 IGB=1,NGB 

CALL YINI (NUT2 ,IH, 1,10) 

NCIGB = I H ( 1 ) 
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NELIC-B = N*MCIGB 

CALL YIN (NUT 2 »V»1 »NELIGB ) 

X 

C SOLUTION FCP. (G) FROM IU* S 'T )*(G) = <B ) . 

DO 27 ICU=I,NC-U 

CALL YINI (NUTU *IH « 1 «iCI 

JSU = IH ( 1 ) 

JFU = I H ( 2 '/ 

NEL1GU = TH ( 3 ) 

CALL YIN (NUTU,V,LCU,LSU+NFLIGU-1) 

DO 37 JB=1*NCIGB 

LBSMI = ( JB-I )*N 

LJJU = LSU-1 

DO 36 JU=JSU t JEU 

LIT JU = LJJU+1 

LJJU = LTTJU+LUV( JU)-1 

IT JU = JU-LUV( JU)+1 

IF ( ITJU .EC. JU) GO TO 36 

LJJUM1 = L J JU— 1 

LGP = LESM1+JU 

LG = LPSM1+ITJU-1 

DO 34 LU=LITJU.LJJUM1 

LG = LG+I 

34 V(LGB) = V(LGB) - V(LU)*V(LG) 

36 CONTINUF 

37 CONTINUF 

SOLUTION FOR (Y) FROM (-D-)*(Y)=(G>. 

LYE = 0 

DO 45 JY=!,NCIGB 
LYS = LYF+1 
LYE = LYS «■ N - 1 
LD = LSD-1 
DO 45 LY=LYS *LYE 
LD = LD + 1 

45 V(LY) = V (LY)/V( LD) 

C 

C SOLUTION FOR (2) FROM (U)*(Z)=(Y). 

C U GROUPS APE OBTAINED IN BACKWARDS ORDER. 
DC 57 IGUX=1 »NGU 
IF (IGUX .EC . 1 ) GO TO 55 
BACKSPACE NUTU 
BACKSPACE NUTU 
CALL YINI (NUTU.IH, 1.10) 

JSU - I H ( 1 ) . 

JEU = IH(2) 

NELIGU = IH ( 3 ) 

CALL YIN (NUTU »V » L SU» LSU+NE LIGU-1 ) 

55 BACKSPACE NUTU 
BACKSPACE NU 1 U 
DO 57 JB=1 T NCIGE 
L2SM1 = (JP-1)*N 
LIT JU = LSU+NELI GU 
DO 56 JUX=JSU t JEU 
JU * JSU+JEU-JUX 
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L JJU - LITJU-1 

LITJU = L J„ -LUV(JU) + 1 

ITJU = JU— LUVC JU ) +1 

IF (ITJU .EO. JU) GO TO 56 

LJJUM1 = LJ JU— 1 

L2 = L2SM1+JU 

L2Y = LZSM1+ IT JU— 1 

DO 54 LU=LITJU,LJJUM1 

LZY = LZY+1 

54 V(LZY) = VIL2Y) - V(LU)*V(L2) 

56 CONTINUE 

57 CONTINUE 

DC 72 1=1,10 
72 IH ( I ) = 0 

IH(1) = NCIGB 

CALL YOUTI (NUT1 ,IH,1,1C) 

CALL YOUT (NUTI ,V,1,NELICB) 

89 CONTINUE 

CONVERT 2 FROM FULL COLUMN NOTATION TO SPARSE NOTATION. 
REWIND NUTI 
REWIND NUT2 

CALL Y INI (NUTI *IH» 1*10) 

IH ( 4 ) - NRP*NCB 

IH( 5 ) = 0 

IH(fc) = 0 

IH(7 ) = 5HWH0LE 

CALL YOUTI (NUT2 ,IH, 1*10) 

JZ = 0 

DO 110 TGB=I f NGP 

CALL YINI (NUTI ,IH, 1,10) 

NC = IH ( 1 ) 

NN2PB = NC*NRE 

CALL YIN (NUTI ,V,l,NNZPB) 

LB = 0 

DO IOC J =1 ,NC 
J2 = JZ+1 
DO 100 T Z = 1 , NR f 
LB = LB+1 

100 LUV(LP) = 1 00000* IZ+JZ 
IH( 1 ) = NNZPP 
IH(2) = LUV(l) 

IH ( 3 ) = LUV(NN2PP) 

CALL YOUTI (N'UIZ ,IH, 1 , 10) 

CALL YOUTI (NU72 ,LUV,I,NN2PB) 

110 CALL YOUT (NUT7 , V,1,NNZPB> 

CALL YLORD (!\M'TZ ,V,LUV,KV,NUT 1 ,NUT2) 

RETURN 

END 
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SUBROUTINE YBTAB (NUTA,NUTB,NUTZ,V,LV,KV,NUT1 ,NUT2) 

DIMENSION V(1),LV(1), MHEAD ( 1C ) 

DATA NIT, NOT/5, 6/ 

C 

C TRIPLE MATRIX PRODUCT FOR SPARSE MATRICES. 

C I B(TRANSPOSE) * A * B = Z ) 

C KV/4 MUST BE ECUAL TO OR GREATER THAN, 

C Cl) NUMBER OF ROWS OF MATRIX B (NRB=NRA=NCA ) 

C AND ( 2 ) NUMBER OF COLUMNS OF MATRIX B. 

C CALLS FORMA SUBROUTINES YIN ,YINI , YLORD , YMULT ,YNOZER,YOUT , 

C YOUTI , YPART , YSYMLH , YSYMUH , ZZBOMB. 

C DEVELOPED BY R A PHILIPPUS. JUNE 1969. 

C LAST REVISION BY RL WOHLEN FOR NASA. MAY 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

C NUTB = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX B IS STORED. 

C NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NERROR EXPLANATION 

C 1 = SIZE LIMITATION EXCEEDED. 

'* 2 = INCOMPATIBLE MATRICES. 

CALL YPART (NUTA ,V, l. V ,K V, NUT1 ) 

CALL YPART (NUTB ,V, LV,KV,NUT1 ) 

C 

C GET (A) HEADER INFORMATION. 

REWIND NUTA 
REWIND NUTB 

CALL YINI (NUTA *MHE AD, 1 , 10 ) 

NR A = MHFAD(l) 

NRASAV = NPA 
NCA = MHE AD ( 2 ) 

NPARTA = MHEAD (3 ) 

NNZ A = MHEAD (4 ) 

IASHAP = MHE AD ( 7 ) 

C GET (B) HEADER INFORMATION. 

CALL YINI (NUTB , MHEAD, 1 ,10) 

NR B = MHEAD ( 1 ) 

NCB = MHE AD( 2) 

NPARTB * MHEAD ( 3 ) 

NNZB = MHEAD (4) 

IBSHAP = MHE AD ( 7 ) 

IADf NS= 100*NNZA/NPA/NC A 
I BDENS = 100*NNZB/NPB/NCB 

IF ( IASHAP. EG. 5HLPWEP ) IADFNS=100*( ?*(NNZ A-NR AKNRA) /NRA/NCA 
IF (IASHAP. EG. SHUPRER ) I ADENS= 100* ( ?* (NNZ A-NR A >+NRA ) /NR A/NCA 
IF ( IBSHAP. EC. 5HL0WFR ) IBD ENS=100* ( 2* (NNZ B -NRB )+NRB ) /NR8/NCB 
IF (IBSHAP. FO. SHUPPEP > I PDENS= ICO* ( ?* (NNZP -NRP )+NFB )/NRB/NCB 
IP (NNZA.EC'.O .OR. NNZB.FO.O) GO TO 95 
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C 

t REWIND NUTZ 

MHEAD(l) - NCB 
MHFAD(2) = N«B 

IF (IBSHAP.EG.5HWH0LE ) MHEAD (5) = 0 
IF I IBSHAP.EO.f^VHOLE ) MHEAD(6) = 0 
CALL YOUTI (NUTZ ,MHE AD, 1,10) 

C 

DO 10 I=1,NPARTB 

CALL YINI (NUTB »MHE AD, 1 , 10 ) 

NNZP = HHEAD(l) 

LFELP = MHEAD( 2 1 
LLELP = MHEAD(3) 

CALL YINI (NUTB ,LV,1, NNZP) 

CALL YIN (NUTB, V, 1 , NNZP) 

IF (IBSHAP.NE.5HWH0LE ) GO TO 6 
C 

DO 5 J=1,NNZP 
K=LVIJ)/1 00000 

5 LV( J)=100000*(LV( J)-100000*K)+K 

6 MHEAD (2 ) = LV(1) 

MHEADC3) = LVCNNZP) 

CALL YOUTI (NUTZ , MHEAD, 1,10) 

CALL YOUTI (NUTZ, LV, 1, NNZP) 

10 CALL YOUT (NUTZ ,V, 1 , NNZP ) 

CALL YMULT (NUTZ ,NUTA »NUT1 »V»LV»KV »NUT2 ) 

IF (IASHAP.E0.5HWH0LE ) GO TO 85 

SYMMETRY OF A IS USED FROM HERE TO STATEMENT 85 
NPARTZ=0 
N NZZ=0 
NRFC=0 

CALL YLORD (NUTB ,V,L V,KV,NUT2 , NUTZ ) 

REWIND NUTl 

CALL YINI (NUTl ,MHE AD , 1 , 10 ) 

NRA = MHEAD ( 1 ) 

NCA = MHEAD l 2) 

NPARTA = MHEAD(3) 

NNZ A = MHEAD (4 ) 

REWIND NUTB 

CALL YINI (NUTB »MHE AD, I , 10 ) 

NRB * MHEAD ( 1 ) 

NCB = MHE AD ( 2 ) 

NPARTB = MHEAD (3) 

NNZB = MHEAD (A) 

ISHAP = MHFAD( 7) 

IF (ISHAP. E0.5HWHOLE .OR. ISHAP. EC.4HDIAG ) GO TO 15 
IF (ISHAP. F0.5HL0WER) CALL YSYMUH ( NUTB ,V , LV, KV, NUT2 , NUTZ ) 
IF (ISHAP. FC.5HUPPFR) CALL Y_ YM LH (NUTB ,V, LV,KV,NUT2, NUTZ ) 
REWIND NUTB 

CALL YINI (NUTf. ,MHE AD, 1 , 10 ) 

NRB = MHE AD (I ) 

NCB = MHE AD (2 ) 
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NPARTB = MHEAD ( 3 ) 

NNZB * MHEAD (4) 

15 IF (NNZA.EQ.C .OR. NNZB.EQ.O) GO TO 70 

IF (NCA.GT.KV/4 .CP. NCB.GT.KV/4) GO TO 999 

IF (NRB.NE.NCA) GO TO 999 
C 

IZ=0 

LPBS=KV/4+l 
LPBF=LPBS-I 
LC$=KV/2+l 
LCE=LCS-1+NCB 
LCC$=LCF+1 
LCCE=LCE 
NNZ=KV— LCCS+l 
REWIND NUT? 

C 

DO 20 I=LCS»LCE 
20 V(I)=0. 

C 

DO 55 1=1 tNPARTA 

CALL YINI (NUT 1 »MHE AD* 1*10) 

NNZPA = MHEAD ( 1 ) 

LFELPA = MHEADC2) 

LLELPA = MHEAD (3 ) 

CALL YINI (NUT1,LV,1, NNZPA) 

CALL YIN (NUT1,V,1, NNZPA) 

K=LPBS 
ITRP L=C 
REWIND NUTP 

CALL YINI (NUTB»MHEADf 1*10) 

NREAD=0 

C 

DO 50 INA=1, NNZPA 
I A=LV ( INA )/100000 
JA=LV(INA)-100000*IA 

IF (TA.EO.IZ .AND. ITRBL.EQ.l) GO TO 50 
ITRBL=0 

IF (TA.EO.IZ) GO TO 30 
REWIND NUTP 

CALL YINI (NUTB »MHE AD* 1 * 10 ) 

MREAD=0 

C 

DO 25 INC = LCS*LC f 
IF (V(INC).EO.O.) GO TO 25 
LCCE=LCCE +1 
V ( LCCE ) =V ( INC ) 

LV( LCCE )= IZZ+INC-KV/2 
V ( INC ) =0. 

1 1 (LCCE.LT.KV) GO TO 25 
CALL YOUTI (NUT 2 *LV*LCCS»LCCE) 

CALL YCUT (NUT 2*V»LCCS*LCCE) 

NREC=NPFC+1 

NNZZ=NNZZ+NNZ 


NERR0R=1 

NFRR0R-2 
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LCCE=LCF 
25 CONTINUE 


IZ = IA 

I2ZelOOOOC*I2 

K=LPBS 

30 IF (K.LE.LPBE .AND. NREAD.GT.O) CO TO AO 
35 IF (NREAD.FQ.NPAPTB) ITP.BL=1 
IF (ITRBL.EQ.l ) GO TO 50 
CALL YINI (NUTB,MHEAD,1,10) 

NNZPB = MHEADC1) 

LFELPB = MHFAD ( 2 ) 

LLFLPP = MHEADC3) 

LPBE=LPBS-1 +NNZPB 

CALL YINI (NU'3,LVtLPBS,LPBE) 

CALL YIN (NUTB »V»LPBS»LPBE) 

NREAD=NREAl>+1 

K=LPBS 

AO DO A5 1MB=K,LPBE 
K = I N P 

IB=LV ( INB )/l OCOOO 
IF (IB.GT.JA) GO TO 50 
IF (TB.LT.JA) GO TO A5 
JBZ-LVC 1NB)-100000*IB 
IF CJPZ.GT.IA) GO TC A5 
INZ=KV/?+ JBZ 

V(INZ)=V(INZ)+V( 1 NA ) *V ( INB ) 

A5 CONTINUE 

GO TU 35 
50 CONTINUE 

55 CONTINUE 

DO 60 I=LCS ,LCF 
IF (VC I). EC. 0.) GO TO 60 
LCCE=LCf E +1 
VCLCCE )=VCI) 

L V C LCC E ) * IZ Z ♦! -K.V/2 

IF CLCCF.LT. KV) GO TO 6C 

CALL YCUTI (NUT?,LV,LCCS,LCCE ) 

CALL YOUT (NUT 2 »V»LCCS,LCCE) 
NREC=NPFC+1 
LCCE=LCE 
NNZZ=NNZZ+NNZ 
60 CONTINUE 

IF (LCCE.EO.LCE) GO TC 7C 
NNZ'LCCE-LCCS+1 
NN Z Z =VN Z 2 +NN 2 
NREC=NRFC+1 

CALL YOUTI (NUT2.LV, LCCS, LCCE ) 

CALL YOUT (NUT2,V,LCCS,LCCE ) 

70 REWIND NUTZ 
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MHBADC1 ) = N«A 

MHFAD ( 2 ) = NCB 

MHFAD (3) = NR EC 

MHFAD(4) = NN2 2 

MHF AD ( 5 ) = 5H0PDEP 

MHFAD(fe) = C 

MHEAD(7 ) = 5HLCWER 

CALL YOUTI (NUT 2, MHE AD* 1,10 

IF (NNZZ.CT.O) GC TC 75 

DO 71 1 = 1 ,10 

vm = o. 

71 MHFAD(I) = C 

CALL YOUTI (NUTZ »MHE AD * 1 * 1C 1 
CALL YOUTI (N rZ ,MHE AD, 1 , 2 ) 

CALL YCUT (NUTZ, V,l,2) 

RETURN 

75 LZE=KV 
REWIND NUT2 

DO 80 T = 1 » NR E C 

IF (I.FO.NPFC) LZF=LCCS-1+NNZ 

NN2P = LZE-LCCS+1 

CALL YINI (NUT?*LV,LCCS,L2E 1 

CALL YIN (NUT 2 »V»LCCS»LZF) 

MHEAD(l) = NNZP 
MHFAD(2 ) ~ LV(LCCS) 

MHEAD(3) = LV(LZE 1 
DO 76 J=4,1C 

76 MHFAD (C. = 0 

CALL YOUTI (NUTZ ,MHEAD, 1,10) 

CALL YOUTI (NUTZ ,LV» LCCS»LZE ) 

80 CALL YOUT (NUTZ ,V, LCCS ,LZE ) 

CALL YPART (NUTZ ,V,LV,KV, NUT 1 ) 

GO TO 90 

85 CALL YMULT (NUT1 ,NUTB.NUTZ,V,LV,KV,HUT2) 
90 RETURN 

95 REWIND NUTZ 
NN2Z=0 

MHEAD(l) = NCB 

MHEAD(2 ) = NCB 

MHEADO) = NNZZ 

MHEAD(4) = NNZZ 

MHFAD (5 ) = 5H0RDER 

MHEADC6) = KV 

MHEAD(7) = 5HWH0LE 

CALL YOUTI (NUTZ ,MHE AD, 1,10) 

DO 96 1 = 1,10 
VII) = 0, 

96 MHFAD ( I ) = 0 

CALL YOUTI (NUTZ *MHE AD*1 , 10 ) 

CALL YOUTI (NUTZ, MHE AD, 1,2) 

CALL YOUT (NUTZ, V,l,2) 

RETURN 
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C 

999 CALL ZZBOME C5HYBTAB 
END 


.NERROR) 
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SUBROUTINE YDCM3A (NUTA ,NUTU ,NUTD*V ,L V,KV.'\*UTI ,NUT2 ) 

DIMENSION vm,L V<1 )* MHFAD(lO) ,MPHEAD(10) 

DATA EPS/I. E-20/ 

C 

C DECOMPOSE SPARSE MATRIX ( A > TO FORM UPPER TRIANGULAR MATRIX WITH ONES 
C ON DIAGONAL (U) AND DIAGONAL MATPTX (-D-) SUCH THAT 
C (A) = (U)**T * (-D- ) * (U). METHOD ATTRIBUTED TO GAUSS. 

C SPECIAL FORM USED FOR FORMA SUBROUTINE YMOD2A. 

C IF THE WHOLE MATRIX (A) IS INPUT* ONLY THE LOWER HALF IS USED. 

C BAND WIDTH (DIAGONAL UP TO TOP NON-ZERO MUST BZ LESS THAN OR EQUAL 

C TO (KV-N)/2* WHERE N IS MATRIX SIZ r (SQUARE). 

C CALLS FORMA SUBROUTINES YIN *YINI » YLCRD *YCUT .YCUTI *YPART * 

C YTR AMS *ZZ BOMB • 

C DEVELOPED BY R L WOHLEN AND R A PHILIPPUS. NOVEMBER 1971. 

C LAST REVISION PY RL WOHLEN FOP NASA. MAY 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 


c 

NUTA 


LOGICAL NUMBER 

OF 

UTILITY 

TAPE 

ON 

WHICH 

MATRIX A IS STORED. 

c 

NUTU 

— 

LOGICAL NUMBER 

OF 

UTILITY 

TAPE 

ON 

WHICH 

U 

IS 

STORED. 

c 

NUTD 


LOGICAL NUMBER 

OF 

UTILITY 

TAPE 

ON 

WHICH 

D 

IS 

STORED. 

c 

V 


VECTOR WORK SPACE 

• 







c 

LV 

=: 

VECTOR WORK SPACE 

• 







c 

KV 

= 

DIMENSION SIZE 

OF 

V*LV IN 

CALLING 

PROGRAM 

• 


c 

NUT1 

=- 

LOGICAL NUMBER 

OF 

UTILITY 

TAPE 

• 





c 

NUT2 


LOGICAL NUMBER 

OF 

UTILITY 

TAPE 

• 






c 

NERROR EXPLANATION 
J 1 = BANDWIDTH LIMITATION EXCEEDED. 

C 2 = SIZE LIMITATION EXCEEDED (LV). 

C 3 = SIZE LIMITATION EXCEEDED (LV). 

C 4 = MATPIX IS SINGULAR. 

C 

C CONVERT FROM SPARSE TO BAND NOTATION. 

KV04=KV/4 
KV02 = KV/2 
KV02P1 = KV02+1 
LAS=KV04+J 
REWIND NtJT A 

CALL YIN I (NUT A *MHEAD» 1 1 1C ) 

NR A = MHEJD(l) 

KVMN = KV-MRA 
KVMN02 = KVMN/2 
I ASHA P = MHE AD ( 7 ) 

NUTS=NUTA 

IF ( IASHAP.FO.SHUPPEP ) CALL YTPANS ( NUTA,NUTU , V* L V,KV,NUT1 ,NUT? ) 
IF ( IASHAP.EQ. SHOPPER ) NUTS=NUTU 
CALL YLORD <NUTS,V,LV*KV,NUTl f NUT2) 

REWIND NUTS 

CALL YINI (NUT S »MHE AC *1*10) 

REWIND VUT1 

I LV = K VO 4 

JLV = KVC^«NRA 

IE ( JLV.LT.KVO?) JLV=KVC? 

JLVS = JLV 
KP = 1 
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LAM AX = LAS-1 +KVMN02 
LAE = LAS 
JS = 1 
NGROUP = 0 
L AS 1 = KV04 
DO 5 I =LA S *KV 
LV( I ) = 0 
5 vm=o. 

NN22 = 0 

NPARTA = MHEADC3 ) 

NROWS = 1 

DO 20 I=1*NPARTA 

CALL YINI (NUTStMPHEAD.l ,10) 

NN2PA = MPHEADd ) 

CALL YINI (NUTS,LV,1,NN2PA) 

CALL YIN (NL'TS.V.l *NN2PA) 

DO 20 J=1,NNZPA 
IA=LV( J 1/100 one 
JA=LVI J )-10CC00*I A 
IF (IA.LT.JA) CO TO 20 
IF (IA.EC.KP) GO TO 15 
LAS! = LAE 
LAE = LAE+IA-JA+1 
NELR = KP-JS+1 

NERR0R=1 

IF (NFLR.CT.KVMN 02) GO TO 999 
NN22 = NNZZ+NELR 
KP = K P+1 
JS = JA 

NROWS = NROWS+1 

ILV = ILV+i 

LV(ILV) = NELR 

IF (LAE.LE.LAMAX) GO TO 15 

JLV = JLV+1 

NERR0R=2 

IF (JLV.GT.KV) GO TO 999 

NROWS = NROWS-1 

LV(JLV) = NROWS 

NROWS = I 

LAE = LAE— IA+JA— 1 

NGP OOP = NGP OU P+T 

CALL YCUT (NUT1 ,V, LAS, LAE) 

DO 10 L=LAS*LAE 
10 V( L ) =0. 

LAS! = KVC4 

LAE = KVC4+IA-JA+! 

KP = IA 

15 LA s LAS1 JA-JS + 1 
VILA )=V( J ) 

20 CONTINUE 

IF (LAS.GT.LAE ) GO TO 65 
NGROUP = NGROUP+1 
ILV = ILV+1 
LV(ILV) = KP-JS+1 


NERR0R=3 



IF (LV(ILV).GT .KVMN02 ) GC TO 999 
NN22 = NN22+LVI IlV) 

JLV = JLV+1 

IF (JLV.GT.KV) CO TO 999 

LV(JLV) - NROWS 

CALL YOUT (NUT1,V,LAS,LAE) 

65 00 30 I=1,NRA 
30 LVfI> = L V( KV04+I ) 

DO 40 1=1, NCR CUP 
40 LVCKV02 + I ) = LV( JIVS*- 1 ) 

C 

C DECOMPOSITION*. 

C D IN VII THRU N). A,U GROUP A START AT V(N*1) 
C A,U GROUP *• START AT V ( N ♦ 1 ♦ ( K V-N )/ 2 ) • 

C LV(I).I=1,N TS NUMBER CF ELEMENTS IN COLUMN I 
C LV(KV/2-*lG) IS NUMBER OF COLUMNS IN GROUP IG. 
N = NRA 
NG = NGRCUP 
LSGA = N+l 

LSGB = LSCA + (KV-N1/2 
REWIND NUTU 
MHEAD(l) = N 
MHEAD(2) = N 
MHEADI3) = NG 

CALL YOUT I (NUTU ,MHE AD ,1,10) 

CALL YOUTI (NUTU ,LV , 1 , N ) 

JEGA = 0 
DO 195 IGA= 1 ,NG 
REWIND NUT1 
REWIND NUT2 
NUTP = NUT1 
NUTQ = NUT? 

IF (2* (IGA/2) .EC. IGA) NUTP=NUT2 
IF (NUTP .FC. NUT2) NUTQ=NUT1 

C OPERATE ON GROUP A ONLY. 

NCGA = LV (KV02+IGA) 

JSGA = JEGA+1 
JEGA = JSGA+NCGA-1 
LEGA = LSGA— 1 
DO 101 J= JSGA, JEGA. 

101 LEGA = LFGA ♦ LV(J) 

CALL YIN (NUTP,V, LSGA, LEGA) 

LJJ = LSGA-1 

DO 140 J= JSGA, JEGA 

JM1 = J-l 

ITOPJ = J— LV( J ) + l 

L1T0PJ = LJJ-H 

LJJ = LITOPJ+LV( J )-l 

IF (J ,F0. JSGA) GO TO 134 

IF (ITOPJ .EC. J) GO TO 13^ 

ISTART = ITOPJ 
L I J = L ITOPJ— 1 

IF (ITOPJ .CF. JSGA) GO TO 105 

ISTART = JSGA 

LI J = LITOP J+JSGA-ITOP J-l 
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105 LITOPI = LSGA 

IF (ISTART .EO. JSGA) &C TO 110 
I SMI =■ ISTART— 1 
DO 107 1= JSGA, ISM1 
107 LITOPI = LITOPI+L V( I ) 

110 DO 12P I=ISTART,JM1 
LIJ = LIJ+1 
IM1 = I-I 
ITOPI = I-LVin + 1 
IF (ITOPI .LT. I TOPJ ) GO TO 115 
KST ART = ITOPI 

IF (I .FO. KSTART) GO TO 125 
LKI - L ITOPI— 1 
LKJ = LITCPJ+ITOPI— ITCPJ— 1 
GO TO 120 

115 KSTART * I TOPJ 

IF (I .FO. KSTART) GO TO 125 
LKI = LITOPI+ITOP J-IT0P1— 1 
LKJ = LITOPJ-1 
120 DO 122 K=KSTART,IM1 
LKI = LKI +1 
LKJ = LKJ+1 

122 V(LIJi = V(LIJ) - V( K )*V( LKI )*VI LKJ ) 

125 V(LIJ) = V(LTJ )/V (I I 
128 LITOPI = LITOPI+LVU) 

134 V(J) = V(LJJ) 

IP (ITOPJ .10. J) GO TO 139 
LKJ = LITOPJ-1 
DO 138 K- ITCPJ , JM 1 
LKJ = LKJtl 

138 V(J> = V(J) - V( K )*V ( LKJ )**2 

139 NERROR=A 
IF (ABS(V(J) ).LT.EPSI GO TO 999 

140 V(LJJ) = 1.0 

C GROUP A OPFRATE ON GROUP B. 

C I COLUMNS ARE IN GROUP A* J COLUMNS IN GROUP B • 

IF (IGA .EO. NG) GO TO 194 

IC-AP1 = IGA+1 

JEGB = JFGA 

DO 192 IGR= 1GAP1 *NG 

NCGB = LV(KVO?-HGB) 

JSGB - JEGB + 1 
JEGB = JS GE+NCCfi-1 
LEGB = LSGB-1 
DO 151 J=JSGB,JEGB 
151 LEGB = LEGP+LVIJ) 

CALI YIN (NUTP ,V»L SGB ,LEGB ) 

LJJ = LSGP-1 

DO 190 J= JSGB, JEGB 

JM 1 = J-l 

ITOPJ r J-L V( J ) + l 

LITOPJ = LJJ + 1 

LJJ = LITCPJ+LV( J )-l 

IF (ITOPJ .GT. J EGA ) GO TO 1«0 

ISTART = ITOPJ 
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LI J = LITOPJ-1 

IF CTTOPJ .GE. JSGA ) GO TO 155 
I START = JSGA 
LI J = L ITCPJ+JSG A-I TOP J— 1 
155 LITOPI = LSGA 

IF (ISTART .EQ. JSGA) GO TO 160 
I SMI = ISTART- 1 
00 157 I=JSGA? ISMl 
157 LITOPI = LITOPI + LVC I ) 

160 DO 178 I=ISTART t J EGA 
LIJ = LIJ+1 
IM1 = 1-1 
ITPPT = I-LV(I)+1 
IF ( ITOPI .LT. I TCP J ) GO TO 165 
KSTART = ITOPI 

IF (I .FO. KSTART) GO TO 175 
LKI = L ITOPI— 1 
LKJ = LIT0PJ+IT0P1-IT0PJ-1 
GO TO 170 
165 KSTART = ITCPJ 

IF (I .EC. KSTART) GO TO 175 
LKI = LITPPI+ITOP J-ITOPI-1 
LKJ = LITOPJ-1 
170 DO 17? K=KSTART,IM1 
LKI = LKI +1 
LKJ = LKJ+1 

172 VCLIJ) = V(LIJ) - VCK )*VCLKI )*VCLKJ) 
175 VCLIJ) = VCLIJ )/Vli ) 

178 LITOPI - LITCPI+L VC I ) 

190 CONTINUF 

192 CALL YOUT CNUTQ ,V , LSGB.LEGB ) 

194 MPHEAD Cl) = JSGA 
MPHEADC2) = JEGA 
MPHEADC3) = LECA-LSGA + 1 

CALL YOUTI (NllTU *MPHE AD*1 *1 0 ) 

195 CALL YOUT CNUTU»V* L$GA,LEGA ) 

REWIND NUTD 

CALL YOUT CNUTD t V,l»N) 

RETURN 

999 CALL ZZBOMP C 6NYDCM2A , NERROR ) 

END 
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SUBROUTINE YDCOM3 (NUTA, NUTU ,NUTO* V,LV,KV , NUT I ,NUT2) 

DIMENSION V(1),LV(1), MHEAD ( 10 ) *MPHE AD ( IP) 

DATA FPS/l.E— 2C/ 

C 

C DECOMPOSE SPARSE MATRIX (A) TO FORM UPPFR TRIANGULAR MATRIX WITH ONES 
C ON DIAGONAL III) AND DIAGONAL MATRIX (-D-) SUCH THAT 
C (A) = (U)**T * (-D-) * (UK METHOD ATTRIBUTED TO GAUSS. 

C IF THE WHOLE MATRIX I A) IS INPUT, ONLY THE LOWER HALF IS USED. 

C BAND WIDTH (DIAGONAL UP TO TOP NON-ZERO) MUST BE LESS THAN OR EQUAL 

C TO ( KV-N )/2» WHERE N IS MATRIX SIZE (SGUARF). 

C CALLS FORMA SUBROUTINES YIN ,YINI ,YLORD ,YCUT ,YOUTI ,YPART , 

C YTRANS,ZZBOMB. 

C DEVELOPED BY R L WOHLEN AND R A PHILIPPUS. DECEMBER 1972. 

C LAST REVISION BY WA BENF1ELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS (AIL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

C NUTU = LOGICAL NUMBER OF UTILITY TAPE ON WHICH U IS STORED. 

C NUTD = LOGICAL NUMBER OF UTILITY TAPE ON WHICH D IS STORED. 

C V = VECTOR WORK SPACE. 

C LV - VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPF. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NERROR EXPLANATION 

F 1 = BANDWIDTH LIMITATION EXCEEDED. 

2 = SIZE LIMITATION FXCEEDED (LV). 

3 = BANDWIDTH LIMITATION EXCEEDED. 

4 = MATRIX IS SINGULAR. 

CONVERT FROM SPARSE TO BAND NOTATION. 

KV04=KV/4 
KV02 = KV/2 
KVP2P1 = KV02+1 
LAS=KV04+1 
REWIND NUT A 

CALL YINI (NUTA,MHEAD,1,10 ) 

NRA = MHFAD(l) 

KVMN = KV-NRA 
KVMN02 = KVMN/2 
I ASHA P - MHE AD ( 7 ) 

NUTS=NUTA 

IF ( IASHAP.F0.5HUPPFR ) CALL YTPANS ( NUTA, NUTU , V, l V,KV ,NUT1 ,NUT2 ) 

IF (IASHAP.FC.5HUFPFR ) NUTS=NUTU 
CALL YLORD (NUTS *V,LV,KV,NUT1 ,NUT2 ) 

REWIND NUTS 

CALL YINI (NUTS , MHEAD, 1,10) 

REWIND NUT! 

ILV = KVC4 
JLV = KV04+NRA 
IF ( JLV.LT.KVO?) JLV-KV02 
JL VS = JLV 
KP - 1 

LAMAX = LAS-1+KVMN02 
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LAE = LAS 
JS = 1 
NGROUP = 0 
L AS 1 = KVOA 
DC 5 I=LAS,KV 
LVI I ) = C 
5 V(I)=C. 

NNZ2 = 0 

NPARTA = MHEADC3) 

NROWS = 1 

DO 20 1=1, NPARTA 

CALL YINI (NUT S »MPH E AD ,1 *10 ) 

NNZPA = MPHEADd I 

CALL YTNI (NUT$,LV,l,Ni'J2PA) 

CALL YIN (NUTS, V,l, NNZPA) 

DO 20 J=l, NNZPA 
IA=LV( J )/10CP0C 
JA=LV( J )— 100000*1 A 
IF (IA.LT.JA) CO TO 20 
IF (IA.EO.KP) GO TO 15 
LAS1 = LAE 
LAE = LAE+IA-JA+1 
NELR = KP-JS+1 

NERR0R=1 

IF (NELF.CT.KVMN02) GO TO 999 
NNZZ = NNZZ+NELR 
KP - KP+1 
JS * JA 

NROWS = NROWS+1 

ILV = ILV+1 

LV(ILV) - NELR 

IF (LAE.LE.LAMAX) GO TO 15 

JLV = JLV+1 

NERR0R=2 

IF (JLV.GT.KV) GO TO 999 

NROWS = NROWS-1 

LV(JLV) = NROWS 

NROWS = 1 

LAE = LAE-IA+JA-1 

NGROUP = NGROUP^l 

CALL YOUT (NUT1 ,V, LAS, LAE) 

DC 1C L=LAS»LAE 
10 V(L)=0. 

LAS1 = KVOA 

LAE = KV0A4 I A-JA+1 

KP = IA 

15 LA = LAS1 + JA-JS+1 
VILA )=V ( J ) 

20 CONTINUE 

IF (LAS. GT. LAE) GO TO 65 
NGROUP = NGROUP 4 1 
ILV = ILV41 
LV(ILV) = KP-JS+1 

NERR0R=3 

IF (LV(ILV).GT .KVMN02 ) GO TO 999 
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NNZZ = NNZZ+LVIILV) 

JLV = JLV+1 

IF IJLV.GT.KV) GC TO 999 

LVULV) = NRCWS 

CALL YCUT (NUT1 »V»LAS*LAE) 

65 DO 30 1 = 1 *NRA 
30 LVCI) = L VI K VC* ■» I ) 

DO 4-0 1=1 ,NGROUP 
40 LVIKV02+I ) = LVIJLVS+I) 

C 

C DECOMPOSITION. 

C D IN VII THRU N). A»U GROUP A START AT VlN+1). 
C A » U GROUP P START AT VIN + l*IKV-N)/2 ) . 

C LVCI), 1=1, N IS NUMBER OF ELEMENTS IN COLUMN I. 
C LVIK V/2+IG ) IS NUMBER OF COLUMNS IN GROUP IG. 

N = NRA 
NG = NGROUP 
LSGA = N4l 

LSGE = LSGA 4 IKV-N>/2 

REWIND NUTD 

JEGA = 0 

DO 195 IGA=1 »NG 

REWIND NUT1 

REWIND NUT? 

NUTP = NUT1 
NUTQ = NUT2 

IF 12*1 IGA/2 ) .EO. IGA) NUTP=NUT2 
IF I NUTP .EG. NUT2) NUTC=NUT1 

C OPERATE ON GROUP A ONLY. 

NCGA = LV IKVC2 + IC- A) 

JSGA = JEGA+I 
JEGA = JSGA4NCGA-1 
LEGA = LSGA— 1 
DC 101 J= JSGA, JEGA 
101 LEGA = LEG:. 4 LVIJ) 

CALL YIN l NUTP ,V, LSGA, LEG A) 

LJJ = LSGA-1 

DO 140 J=JSGA» JEGA 

JMI - J-I 

ITOPJ = J-LVIJ)4l 

LITOPJ = LJJ+1 

LJJ = LITCPJ4LVI J)-l 

IF IJ .EG. JSGA) GO TO 134 

IF ( I TOP J .EO. J) GD TO 134 

ISTART = ITOPJ 

LI J = LTTOPJ-1 

IF ( ITOPJ .GE. JSGA) GO TC 105 
ISTART = JSGA 
LIJ = L ITOP J+JSGA-1T0P J-l 
105 LITOPI = LSGA 

IF IISTART .EO. JSGA) GO TC 110 
ISM1 = TSTART-I 
DO 107 1 = JSGA, ISM 1 
107 LITOPI = IITOPULVII) 

110 DO 128 I=ISTART,JM1 
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tIJ * LIJ+1 
S - V(LIJ) 

IM1 = 1-1 

ITOPI = I— LV ( I )+l 

IF (ITOPI .LT. I TOP J ) GO TO 115 

KSTART = ITOPI 

IF (I .FQ. KSTART) GO TP 125 
LKI = LITOPI-1 
LKJ * LITOPJ+ITOPI-ITOPJ-1 
GO TO 120 

115 KSTART = 1TPPJ 

IF (I .EO. KSTART) GO TO 125 
LKI =: L ITOPI+ITOP J-IT0P1— 1 
LKJ = L ITOP J— 1 
120 DO 122 K=KSTAPTtIMl 
LKI * IKI+1 
LKJ =: LKJ+1 

122 S = S-VCK )*V(LKI )*V(LKJ) 

125 V(LIJ) = S/V(I) 

128 LITOPI = LITOPI+LV( I ) 

134 V(J> = V(LJJ) 

IF ( ITOP J .EO. J) GO TO 139 
LKJ = LITCPJ-1 
DO 138 K=IT0PJ,JM1 
LKJ = LKJ+1 

138 V(J) * V(J) - V(K)*V(LKJ)**2 

139 NERR0R=4 
IF (ABS(V(J) ).LT.EPS) GO TO 999 

140 V(LJJ) = 1 .0 

C GROUP A OPFRATE ON GROUP B. 

C I COLUMNS ARE IN GPOUP A t J COLUMNS IN GROUP B. 

IF (IGA .EO. NG) GO TO 195 

IGAPl = IGA + 1 

JEGB = JFGA 

DO 192 IGB= IGAPl ,NG 

NCGP = LV(KV02+I GB ) 

JSGB = JEGB+1 
JEGB = JSGB+NCGB-1 
LEGB = L5GB-1 
Dn 151 J= JSGE * JFGB 
151 LFGB = LFGF+LV(J) 

CALL YIN (NUTP f V»L SGB » LFGB ) 

LJJ = LSGB-1 

DO 190 J= JSGB* JEGB 

JM1 = J-l 

ITOP J = J— LV(J)+1 

L1T0PJ = LJJ+1 

LJJ = l TTOP J+LV ( J )-l 

IF (ITOFJ .GT. JFGA) GO TO 190 

ISTART = ITCPJ 

LI J = LITOPJ-I 

IF (ITOPJ «GF. JSGA) GO TO 155 
ISTART = JSGA 
LTJ = LITPPJ+JSCA-ITOPJ-1 
155 LITOPI = LSGA 



IF (ISTAPT .EC. JSCA) GO TO 16C 
ISM1 = ISTART-! 

DO 157 1= JSGAt ISM1 
157 L1TOPI = L ITOP I + LV( I ) 

160 DO 178 I=ISTART,JEGA 
LI J = LIJ+1 
S = V(LIJ < 

IH1 = 1-1 

ITOPI = I-Lvm + 1 

IF (ITOPI «LT» ITOPJ ) GO TO 165 

KSTART = ITOPI 

IF CT .EC. KSTART) GO TO 175 
LKI = L ITOPI— 1 

LKJ = LITCPJ+ITOPI— ITOPJ-1 

GO TO 17C 

165 KSTART = ITOPJ 

IF (I .EC. KSTART) GO TO 175 
LKI = L ITOPI+ITO PJ— ITOPI— 1 
LKJ = LITCPJ-1 
170 DO 172 K=KSTART,IM1 
LKI = LKI+I 

LKJ = LKJ+1 

172 S = S— V ( K ) *V (LKI ) *V ( L K J ) 

175 V(LIJ) = S/V(I ) 

178 LITOPI = LITOPI+LVCI) 

190 CONTINUE 

192 CALL YOUT (NUTQ ,V, L SGB , LEGB ) 

195 CALL YOUT (NUTD ,V,L SGA.LEGA ) 

C 

C CONVERT FROM PaND TO SPARSE NOTATION. 
REWIND NUT2 

CALL YOUT (NUT 2 tV» 1 »N ) 

REWIND NUTU 
REWIND NU T D 
LVGS = KV-NGRCiUP 
LVR * LVGS 

DO 202 IGR0UP=1 » NGP.OUP 
LVR = LVR +1 

202 LV(LVR) = LV(KV02+IGR0UP ) 

LS = LVGS— N 
LVE = LS 
DO 2C4 1=1 ,N 
LVE = LVE +1 
204 LV(LVF) - LV ( I ) 

KVMAX = KV/4 

IF (KVMAX. GT.LS) KVMAX=LS 

MHEAD ( 1 ) = N 

MHEAD( 2 ) = N 

MHEAD ( 3 ) = NGRO'JP 

MHEAD (4 ) = NN22 

MHEAD ( 5 ) = 0 

MHE AD ( 6 ) = 0 

MHE AD ( 7 ) = 5HWH0LE 

CALL YOUTI (NUTU »MHE AD* 1 ♦ 1C ) 

LVI = 0 
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LVR = LVGS 
LVF = LS 
LVEP = LS 
IZ = 0 

DO 250 IGROUP=l »NGROUP 
LVR = LVR +1 
LZ = 0 

NROWS = LV(LVR) 

NELG = 0 

DO 206 1=1, NROWS 

LVE = LVE+1 

2 06 NELG = NELG+LV(LVE) 

CALL YIN L’JUTD ,V» 1 ,NELG ) 

DO 2 OP 1=1, NROWS 
1Z = IZ+1 
LVEP = LVEP+I 
JS = I Z-LV( LVEP) +1 
DO 208 J2=JS,IZ 
LZ = LZ+1 

208 LV(LZ) = 10C000*JZ+IZ 
MPHEAD(l) = LZ 
MPHFAD ( 2 ) = LVll) 

M PH F AD ( 3 ) = LV(LZ) 

CALL YOUTI (NUTU ,MPHE AD ,1,10) 
CALL YOUTI (NUTU,LV, 1 ,LZ ) 

CALL YOUT (NUTU ,V , 1 , LZ ) 

250 CONTINUE 

REWIND NUT? 

REWIND NUTD 

CALL YIN (NUT2 »V» 1 »N ) 

DO 260 1=1, N 

260 LVIII = 100000*1+1 
MHE AD ( 3 ) = 1 
MHFAD(A) = N 
MHEAD (7) = 4HDIAG 
CALL YOUTI (NUTD »MHE AD, 1 , 10 ) 
MPHEAD ( 1 ) = N 
MPHEADf 2 ) = LV( I ) 

MPHEAD ( ? 1 = LV(N ) 

CALL Y0l-.11 (NUTD, MPHEAD, 1,10) 
CALL YOUTI (NUTD ,LV, 1 ,N ) 

CALL YOUT (NUTD ,V » 1 , N ) 

CALL YPART (NUTU ,V, L V ,KV *NUT1 ) 
RETURN 

999 CALL ZZB0M6 (6HYDCOM3 ,NFRROR ) 
END 
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SUBROUTINE YDIAG (NUT A, NUTZ , V , l V,KV ) 

DIMENSION VI! ) ,L V (1 ) , MHE AD (1 0 ) 

DATA NIT tNOT/5 ,6/ 

DIAGONALIZE A SPARSE VECTOR (ROW OP COLUMN) INTO A SQUARE MATRIX. 
CALLS FORMA SUBROUTI‘€S YIN ,VINI ,YOUT ,YOU?l ,2ZBOMB. 
DEVFLCPFD BY W A BENFIELD. OCTOBER 1970. 

LAST RFVISION BY WA BENFIELD. MARCH 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX 2 IS S r OREO« 

V = VECTOR WORK SPACE. 

LV = VECTOR WC'£ K SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

NERROP EXPLANATION 

1 = MATRIX IS NOT A VECTOR. 

2 = SIZE LIMITATION EXCEEDED. 

REWIND NUTZ 
REWIND NU T A 

CALL YINI ( NUT A, MHE AD ,1 ,10) 

NR A = MHFADU) 

NC A = MHE AD ( ? ) 

NP ART = MHE AD| 3) 

MHEAD (7 ) = 4 HD I AG 

n R OR — 1 

IF (NRA.NF.I .AND. NCA.NE.I) GO TO 999 

IF (N D A .EC. 1 ) NRA = NC A 

IF (NCA .EC. 1) NCA = NRA 

MHEAD ( 1 ) NRA 

MHEAO (2) = NCA 

CALL YOUTI (NUTZ .MHEAD. 1 » 10 ) 

NERR OR =2 

DO 30 1 = 1 , N P AR 7 

CALL YINI (NUTA tMriEAD.lt 10) 

NN2P = MH E AD (1 ) 

IF (NNZP ,GT . KV) GO TO 999 
CALL YINI (NUTAtLV.l tNNZP) 

CALL YIN (NUTA. V,l, NNZP) 

DC 20 K=1 .NNZP 
I A * LV(K)/1C0000 
IF (IA.F0.1 ) GO TO 10 
LV(K ) = 1 00000*1 A ♦ I A 
GO TO PC 

10 JA = LV(K ) —100000*1 A 

LV(K) = 1 OCOOC^j A ♦ JA 
20 CCNTINUF 

MHFADm = LV( 1 ) 

MHE AD ( 3 ) = LV(NNZP) 

CALL YOUTI (NUTZ .MHEAD, 1 .10 ) 

CALL YOUTI (NUTZ, LV.l .NNZP) 

30 CALL YOUT (NUT Z ,V , 1 , NNZP ) 



YDIAG — 2/ 2 


C 

RETURN 

999 CALL ZZBCMB <5HYD!AC 
END 


tNEFRC R ) 
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SUBROUTINE YDISA (NUT A , I RA* JC A ,NUTZ ,NRZ ,NC 2 , V, LV , KV,NUT1 ) 
DIMENSION V(1),LV(1), MHFAD ( 1C ) 

DATA NlTtNOT/5,6/ 

C 

C SPARSE MATRIX DISASSEMBLY* (MATRIX Z FROM MATRIX A), 

C CALLS FORMA SUBROUTINES YIN ,YINI » YLORD ,YCUT ,YCUTI * YPART * 

C ZZBOMB. 

C DEVELOPED BY R A PHILIPPUS. FEBRUARY 1«>7C. 

C LAST RFVISICN BY WA EENF1ELD. MARCH l*J7fc. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

C IRA = ROW NUKPEP IN MATRIX A OF FIRST ROW OF MATRIX Z. 

C JCA = COLUMN NUMEFP IN MATRIX A OF FTRST COLUMN OF MATRIX Z. 

C NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

C NRZ = NUMB EP OF ROWS IN MATRIX Z. 

C NCZ = NUMPEP CF COLUMNS IN MATRIX Z. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NERROP EXPLANATION 

C 1 = LOOKING FOR DATA OUTSIDE OR MATRIX A. 

C 

CALL YLCPD (NUTA ,V, LV ,KV,NUT1 , NUTZ ) 

REWIND NUTA 

CALL YINI (NUTA ,MHEAD, 1,10) 

NR A = MHFAD(l) 

NC A = MFE AD ( 2 ) 

NPARTA = M HR AO ( 3 ) 

NNZA = MhFAD(M 
ISHAP = MHF AD ( 7 ) 

IF (IPA.EC.JCA .OR. ISHAP. E0.5HWH0LF) GO TO 5 

IF (ISHAP. FO.SHLOWER) CALL YSYMUH ( NUTA ,V , LV, KV,N'JT1 , NUTZ ) 

IF (ISHAP. EC. 5HUPPEP) CALL YSYMLH ( NUTA ,V, LV, KV.NUT1 , NUTZ ) 

REWIND NUTA 

CALL YINI (NUT A *MHE AD , 1 , 1C ) 

NR A = MHFAD(I) 

NCA = MHEA0( ?.) 

NPAPTA = MHEAD (“ ) 

NNZA = MHFAD(A) 

ISHAP = F Lr A0<7) 

5 LZ S=KV/A+ 1 
LZ F=LZS-1 +KV/A 
LZ=IZS-I 
IP AZ = IP A- 1+NF-Z 
JCAZ= JCA— 1*NCZ 

NERROR=l 

IF (IPAZ.GT.NRA .fR. JC A2 .GT .NCA ) GO TO 999 

I J7«l OCCCO* ( TRA-1 )*JC A— 1 

NNZZ=C 

N P A P T Z = C 

PEWIND NUTT 


C 
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DO 30 1=1 ,NPARTA 

CALL YINI (NL'TA,MHEAD,1 ,10) 

NNZPA = MFFAD(l) 

LFFLPA = MHEAD(2) 

LLELPA = MHEADO 1 

CALL YINI (NUTA,LV,1, NNZPA) 

CALL YIN (NUTA,V,1 , NNZPA) 

I AF=LFFLP A/100000 

IF (IAF.GT.IRA2 .AND. I.LT.NPARTA) CO TO 3C 
IAL=LLELPA/10CC00 
IF (IAL.LT.IRA) GO TO 30 
C 

DO 20 J=1 tNNZPA 
IA=LV( J)/ ICC 000 

IF (IA.LT.IRA .OR. IA.GT.IRAZ) GO TO 10 
JA=LV( J)-1C0C0C*IA 

IF (JA.LT.JCA .OP. JA.GT.JCA2) GO TO 10 
IF (V(J).FG.O.) CC TO 13 
L2=LZ*1 
V(LZ )=V( J ) 

LV(LZ)=LV(J)-IJZ 

NNZZ=NNZZ*I 

10 IF (L2.GF.LZE) GO TO 15 

IF (I.LT.NPARTA -OR. J.LT.NN2PA) GO TO 20 
15 N=LZ— LZ5+1 

IF (L2.LT.L2S) GO TO 2C 
MHFAD(l) = N 
MHEADC2) = LV(LZT) 

MHE AD (3 ) = LV(LZ) 

CALL YPUTI (NUT1 ,MHEAD,1 ,1C ) 

CALL YOUTI (Nt!Tl,LV,L2S*L2) 

CALL YOL’T (F'UTltVtLZSfLZ) 

LZ=LZS-1 
NPART2=NPART2+1 
20 CONTINUE 
30 CONTINUE 
C 

RFW1N0 NUT1 
REWIND NUT 2 
MHE AD (1) = NR2 
MHE AD ( 2 ) = NCZ 
MHFAD ( 3 ) = NPARTZ 
MHE AD ( A ) = NNZZ 
MHF AD ( 5 ) = 5H0RDFR 
MHE AD ( 7 ) = ISHAP 
CALL YOUTI (NUT ZtMHEAD, 1*10 
C 

DO AO I=1,NPAPTZ 

CALL YINI (NUTl.LV, 1,10) 

CALL YOUTI (NUTZ.LV.l ,10) 

N=LV ( 1 ) 

CALL YINI (NUT 1 ,LV* 1,N ) 

CALL YIN (NUT1,V,1,N) 

CALL YOUTI (NUT 2,LV,1,N) 

AO CALL YCUT (NUTZ,V,I,N) 
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C 

CALI YRaRT (NUT2 ,V,LV,KV,WJT1 ) 
RETURN 
C 

999 CALL 22BCWB (5HYDISA * NERRCR ) 
ENO 
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SUBROUTINE YDTCS (A,NUTA ,NRA ,NCA,KR A ,KC A, V ,LV ,KV , NUT1 ) 

DIMENSION V(II* LVMli AIKPA.I), MHEAO(IO), MPHEAD(IO) 

DATA NIT*NCT/5*fc/ 

CONVERT A MATRIX FROM DENSE NOTATION TO SPARSF NOTATION. 

CALLS FORMA SUBROUTINES YIN *YINI *YOUT ,YOUTI ,YPART .ZZBOM6. 
DEVELOPED BY » A RHILIPPUS. JANUARY 1969. 

LAST REVISION BY WA BENFIELD FOR NASA. MAY 197 fe. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

A = DENSE MATRIX. SIZE ( NRA , NCA ) • 

NUTA = LOGICAL NUMBER OF UTILITY T A p p ON WHICH SPARSE MATRIX A WILL 
FF STCRFD. 

NRA = NUMBER OF ROWS IN A. 

NCA = NUMPER OF COLUMNS IN A. 

KRA = ROW DIMENSION OF A IN CALLING PROGRAM. 

KCA = COLUMN DIMENSION OF A IN CALLING PPCGRAM • 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

N’JTl = LOGICAL NUMPER OF UTILITY TAPE. 

NEPROR EXPLANATION 
1 = SIZE LIMITATION EXCEEDED IKA). 


NERROR=l 

IF (NRA. GT. KRA .CP. NCA. GT. KCA) GO TO 999 
L=0 

REWIND NUT1 
NNZ A=0 
NREC= C 
DO 2 1 = 1 * IC 
MHEAD(I) = C 

2 MPHEAD(I) = 0 

STORE PARTITIONS ON NUT1 TEMPORARILY. 

DO 20 1 = 1* NR A 
DO 1C J=1,NCA 
IE (A(I,J).EO.O. ) GO TO 3 
L=L ♦! 

V(L)=A( I* J) 

LV(L)=100C00*I+J 

3 IF (L.GE.KV/A) GO TO 5 
IF (I.LT.NRA) GO TO 10 
IF (J.LT.NCA) GO TO 1C 
IF (L.FO.O) GO TO 1C 

5 MPHEAD(l) = L 
MPHEAD ( 2 ) - LV ( 1 ) 

MPHEADO) = LVCL) 

CALL YOU T I (NUT) .MPHEAD, 1 *1C ) 

CALL YOUT 1 (NUT! ,LV, 1 ,L ) 

CALL YOUT (NUT 1 *V * 1 * L ) 

NREC = NFFC+1 
NN2 A=NNZA +L 
L = C 
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10 CONTINUF 
20 CONTINUE 

REWIND NUTA 

MHFADU) = NR A 

MHF AD ( 2 ) = NCA 

MHFADO) = NREC 

MHF ADC 4 ) = NNZA 

MHFAD (51 = 5HOPDEP 

MHFADm = 5HWHOLF 

CALL YOUTI (NUTA,MHEAD,1,10) 

REWIND NUT1 

TRANSFER PARTITIONS FROM NUT1 TO NUTA. 
DO 65 I = 1,NFFC 

CALL YINI INUT1,MPHEAD,1,10) 

CALL YOUTI (NUT A »M PH FAD *1 .10 ) 

NNZP = M PHF AD ( 1 ) 

CALL YINI (NUT1 ,LV f 1,NNZP) 

CALL YIN (NUT1 ,V,1 ,NNZP) 

CALL YOUTI (NUTA,LV,1 ,NNZPI 

CALL YOUT (NUTA,V*1,NNZP) 

65 CONTINUE 

CALL YPAPT (NUTA,V»LV,KV*NUTI ) 
PFTURN 

999 CALL Z?eOMB (5HYDYQS ,NERROR| 

END 
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SUBROUTINE YEQUALCNUT A,NUTB tVtLV.KV) 

DIMENSION V<n.LVm,MH(10l 
C 

C THIS SUBROUTINE COPIES A SPARSE MATRIX ON UTILITY FILE NUTA 
C TO UTILITY FILE NUTB. 

C SUBROUTINE ARGUMENTS 

C NUTA - INPUT UTILITY FILE CONTAINING A SPARSE MATRIX 
C NUTB - INPUT UTILITY FILE THE SPARSE MATRIX ON NUTA IS TO BE COPIED 
C ONTO. 

C V INPUT WORK SPACE. 

C LV - INPUT WORK SPACE. 

C KV - INPUT DIMENSION OF V AND LV IN THE CALLING PROGRAM. 

C 

C FORMA SUBROUTINES YIN t YINI t YOUT,YOUTI AND ZZEOMB ARE CALLED. 

C CODED BY JOHN ADM IE F *NASA* DEC FMBFR 197?. 

C 

REWIND NUTA 
REWIND NUTB 

CALL V INI (NUTA *MH»1 » 1 C ) 

CALL Y0UTI(NUTB,MH,1,10) 

NPART=MH( 3) 

NNZ A=MH (4 ) 

DO 1C L=1 tNPART 

CALL YINI (NUTA,MH,1 ,10 ) 

CALL YOUTI(NUTF,MH,l,ir>) 

NNZP-MH ( 1 ) 

NERR0R=1 

IF (NNZP .GT. KV) GO TO 999 
CALL YINI (NUTAtL V »1 »NNZP ) 

CALL YOUTI (NUTP»LV»1» NNZP ) 

CALL YlN(NUTA,V t 1»NNZP) 

CALL YOUT (NUTB »V » 1*NNZP) 

10 CONTINUE 
20 RETURN 

999 CALL ZZBOMP fGHYE GUAL t NERROR ) 

END 
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SUBROUTINE YIN (NUT ,A ,NS,NE ) 

DIMENSION A(l) 

X 

C READ IN BINARY DATA FROM PERIPHERAL UNIT NUT INTO CORE ARZA A. 

C CALLS FORMA SUBROUTINE ZZBOMB. 

C DEVELOPED BY R A PHILIPPUS. NOVEMBER 1971. 

C LAST REVISION BY WA BENFIELD FOP NASA. MAY 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUT = _OGICAL NUMBER OF UTILITY TAPE. 

C A = VECTOR TO EE READ. 

C NS = START LOCATION IN VECTOR A. 

C NE = END LOCATION IN VECTOR A. 

C 

C NERRCR EXPLANATION 

C 1 = START LOCATION GREATER THAN END LOCATION. 

C 2 = END OF FILE ENCOUNTERED. 

C 3 = READ PARITY ERROR. 

C 

NERROR=l 

IF (NS -LE. 0 .OR. NS .GT. NE ) RETURN 

NE RRCR=2 

READ (NUT,ERR=998 t END=999) ( A ( I ) , I=NS ,NE ) 

RETURN 

C 

998 NERR0R=3 

999 CALL ZZBOMB (3HYIN 
END 


.NERRCR) 
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SUBROUTINE YINI (NUT, IA,NS,NE ) 

DIMENSION IA(1) 

READ IN BINARY DATA FROM PERIPHERAL UNIT NUT INTO CORE AREA IA. 
CALLS FO&MA SUBROUTINE ZZBOMB. 

DEVELOPED BY R A PHILIPPUS. MAY 1973. 

LAST REVISION BY WA BENFIELD FOR. NASA. MAY 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUT = LOGICAL NUMBER OF UTILITY TAPE. 

IA = VECTOR TO BE READ . 

NS = START LOCATION IN VECTOR IA. 

NE = END LOCATION IN VECTOR IA. 

NERROR EXPLANATION 

1 = START LOCATION GREATER THAN END LOCATION. 

2 = END OF FILE ENCOUNTERED. 

3 = READ PARITY ERROR. 

NFRRCR=1 

IF (NS .LE. C .OR. NS .GT. NE) RETURN 

NERRCR=2 

READ (NUT, ERR =998 , END =999) ( I A ( I ) , I=NS ,NE ) 

RETURN 
C 

998 

999 CALL ZZBOMB (4HYINI , NERROR ) 

END 


NERROR=3 
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SUBROUTINE YLCPD tNUTA,V f LV*KV,NUTl ,NUT?) 

DIMENSION V(l) ,LV(1),IU(16),ILU6>,MHEAD( 1C),MPHFADU0> ,M2HEAD<10) 
DATA NIT»N0T/5 f 6/ 

C 

C ARRANGE ELEMENT LOCATIONS OF MATRIX A INTO INCREASING ORDER. 

C ARRANGE f LFMENTE OF MATRIX A ACCORDINGLY. 

C CALLS FORMA SUBROUTINES YIN ,YINI ,YOUT ,YOUTI * YPART ,ZZ60MB. 

C DEVELCPFD BY p A PHIL1PPUS. DECFMPER 1«>68. 

C LAST REVISION BY WA BENFIELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

C V - VFCTCP WORK SPACE. 

C LV = VECTOP WORK SPACE, 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NERROR EXPLANATION 

C 1 = TWO LIKE LOCATION NUMBERS ENCOUNTERED. 

C 2 = TWO LIKE LOCATION NUMBERS ENCOUNTERED. 

C 

3001 FORMA! (1H1) 

3002 FORMAT ( A ( I 12 ,?X E 17. 8 ) ) 

C 

CALL YPART (NUT A ,V , L V ,KV,NUT 1 ) 

REWIND NUTA 

CALL YINI (NUTA,MHEAD,1,1C) 

NP ART = MHEAD(2) 

NNZA = MHEAD(A) 

I FORD = MHFAD(E) 

IF ( IFORL'.EO. SHORTER .OR. NNZA .LT .2 ) RETURN 

NR FC=0 

NNZPT=0 

kvoa=kv/a 

NTI =NUT ! 

REWIND NTI 
NDO= (NPA? T— 1 )/4 + l 
C 

DO 225 I JK=1 *NDO 

KJI-0 

LPS = 1 

5 CALL YINI (NUTAtMPHEAD,l »1C ) 

NNZP = MPHFAD(l) 

NNZPT = NNZ P7 +NNZ P 
LPF-LPS— 1 t-NNZP 

CALL YINI (NUTAt! V*LPS,LPE ) 

CALL YIN (NUTA,V,LPS,LPF ) 

KJI=KJ1+1 

LPS=LPE+1 

IF (NNZPT .F0.NN2A ) GO 10 10 
IF (KJT.LT. A) GO TO 5 
C 

C SINGLEION METHOD 
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10 M*1 

LAEM1=IPE— 1 

1=1 

J=LPF 

20 IFU.GE.J) GO TO 170 
110 K=I 

I J= < J+I )/2 

IT=LV( I J) 

IMLVm.LE.IT) GOTO 120 
LV(IJ)=LVCI) 

LVU) = IT 
IT=LV<IJ| 

TC=V<IJ) 

V<IJ1=VII ) 

V( I l=TG 
120 L=J 

IF( LVf J 1 »GE • IT) GO TO 140 
LV(IJ»=LV<J) 

LV(J)=IT 

IT=LVIIJ) 

TG=V(IJ) 

V<IJ)=V(J) 

V(J)=TG 

IF(LVm.LF.IT) GO TO 140 

LVCIJ)=LV(I) 

tvm = lT 

IT=LVI I J) 

TG=VUJ) 

v ( u )= vn ) 

Vf I)=TG 
GO TO 140 
130 LVai=LV(K) 

LV(K)=1TT 

TG=VCL1 

V ( L )= V ( K ) 

V(K)=TG 
140 L = L-l 

IF(LV(L l.GT.ITI GO TO 140 
ITT=IV ( l ) 

150 K=K*1 

IFdVlKl.LT.IT) GO TO 150 
IF(K.LF.L) GO TO 130 
IF(L-I.LE.J-K) GO TO 160 
It ( M ) = I 
IU(M)=L 
I = K 
M = K +1 
GO TO IPO 
160 UIHI-K 
I U( M ) = J 
J=L 
M = M + 1 
GO TO 180 
170 M=K~1 
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IFIM.EO.O) GO TO 210 
I=IL(M) 

J=IU(M ) 

180 IFCJ-I.GF.ll) GO TO 110 
IF (I.EO.l) GO TO 20 
1 = 1-1 
190 1=1+1 

IFtX.FO.J) GO TO 170 
IT=LV( I +1 ) 

IF(LVm.LE.IT) GO TO 190 
TG=VC 1+1 ) 

K=I 

200 LV ( K+l )=LV(K) 

V CK + 1 )=V(K) 

K=K-1 

IF ( IT. LT. LV ( K ) ) GO TO 200 
LV(K+1 ) = I T 

V ( K + l ) = TG 
GO TO 1 90 

210 NFRRCR=1 

DO 215 1=1 tLAEMI 
IF avm.EQ.LVU+m GO TO 900 
215 CONTINUE 

IF (LPE.LE.KV04*3 ) KJ 1= ( LPE-1 l/K V04+1 
L PS = 1 

IF (NPART.fT.4) CO TO 21F 
NPART=(LPF-1 )/i<V04+l 
RFWIND NUT A 
MHEAD(3) = NPART 
MHE AD ( 5 ) = 5HP P DFP 
MKFAD (6 ) = 0 

CALL Y0UT1 (NUTA,MHE AD.ItlO ) 

DC 217 T= 1 f NP A p T 

L p PE=LPS-l+KVC4 

IF (LPPE.GT.LPE) LPPE=LPE 

NNZ=LPPF-LPS+1 

MPHFAD ( 1 I = NNZ 

MPUE AD ( 2 ) = LVtLPS) 

NPHEADO) = LV(LPPE) 

CALL YOUTI (NUTA.MPHE AD , 1 ,10 ) 

CALL YPL'TI (NUTA »LV, LPS,LPPE ) 

CALL YOUT (NUTAfV»LPS»LPPF ) 

217 LP5=LPPE+1 
GO TO 310 

21 P DO 220 J=1,KJI 

IF (LPS.GT.LP? > GO TO 225 
LPPE=LPS-1 +KV04 
IF (LPPF.GT.LPF) LPPF=LPE 
NNZ = LPPF-1 PS 
MPHFAD(I) = NMZ 
MPHFADt?) = LV(LPS) 

MPHEAD(3) = LV(LPPF) 

CALL YOUTI ( N T 1 , MPH EAD»1«1C) 
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CALL YOUTI (NTi ,LV,LPS,LPPF ) 

CALL YOUT (NTl , V,LP S ,LPP E ) 

NREC=NREC+1 
220 LPS=LPPE+1 

225 CONTINUE 

NOW THERE ARE NPEC ORDERED PARTITIONS WRITTEN ON NT1 
REWIND NUTA 
MHFAD ( 3 ) = NR EC 
MHF AD ( 5 ) = 5H0RDFP 
MHEADIfc ) = 0 

CALL YOUTI JNUTA,MHEAD,1,1C) 

NT2=NUT2 


MESHING OPER ATICN 
230 REWIND NT 1 
REWIND NT 2 

CALL YINI (NT1 , MPHE AD , 1 , 10 ) 

NN2P3 = MPHEADd ) 

CALL YINI (NTI t LV t l,NN2Pl) 

CALL YIN (NT1 » V* 1 » NNZP1 ) 

IF (NRFC.EC.l) GO TO 305 

LP2S=NN2P7+1 

NPEC2=0 

DO 300 1=2 1 NREC 

CALL YINI (NTl * M2HE AD* 1 » 10 ) 

NNZP2 = M2HE AD ( 1 ) 

LP?E = LP2S-1 +NNZP 2 

CALL YINI (NT1 ,LV,LP2S,LP2F ) 

CALL YIN (NTl ,V»LP?S tLP?E > 

IF ( L V ( LP2S > .GT. LV(NNZP1>> GO TO ?<>5 

MESH TWO PARTITIONS 
Il=i 

1 2 = NNZ P 1 + 1 
IV=2*KVD4 
I2 = C 

250 IW=IW+1 

NE PR0R=2 

IF (LVtll )-LVUZ ) ) ?65,90R,?55 
255 V(IW)=V(I?I 
LV»IW) = LVU2) 

12' 12 + 1 

IE (12.GT.LP2E) GOTO 275 
GO TO 250 
265 V(IW)=V(I1 ) 

LV(IWI=LV(I1) 

11=11+1 

IF (II .(.T.NNZP1 ) GO TO 285 
fC TO 250 

275 NrLTH=NNZP1-Il+l 
K = L P2F 
L=NNZP1 



o o 


DO 280 J= 1 ,NELTM 
V ( K ) ~V ( L ) 

LV(K)=LV(L) 

K = K - 1 
28C L=L-1 
C 

285 IF (IW.t'C.2*KV04) GC TO 295 
J1=2*KVD4 + 1 
C 

to ?c o J=J1,IW 

IZ=IZ-H 

V(IZ)=V(J) 

290 LV(IZ)-LV(J> 

C 

295 NRFC2=NPFC2+1 

M2HEAD ( 2 1 = LV(LP?S) 

M2KE AD (3) = LVO.P2E) 

CALL VOUTI (NT2,M2HEAD,1»1C) 
CALL YOUTI (NT?,LV,LP2S,LP2E ) 
CALL YOl'T (NT2,V,LP2S,LP2E ) 
300 CONTIWIIF 


C 


c 


ALI NRFC PARTITIONS HAVE PEEN RJ AD FRCP 1 NT1 
MPHEAD(2) = LV(1 ) 

NPHFAD ( 3 ) = LVINN2P1) 

CALL YOl'T I (NUTA ,MPHEAD, 1 , 10 ) 

CALL YOUTI (NUTA,LV» 1 ,NNZP1 ) 

CALL YOL'T (NL TA ,V, 1 ,NNZP1 ) 

NREC=NRFC2 
NTS=NT1 
NT! =NT2 
NT2=NTS 
GC TO 23C 


305 CALL YOUTI 
CALL YOUTI 
CALL YDUT 
310 CALL YPART 
RFTURN 


(NUTA,MPHEAD,1 ,10) 
(NUT A ,LV 1 1 ,NN2P1 ) 
(NUT A ,V» i , NNZP 1 ) 
(NUT A ,V, L V? KV,NUT1 ) 


900 WRITF ( NOT , 30C 1 ) 

WPITF (NOT, 3002) LV ( I ) 

WRITF (NOT, 300?) (L V( 1 1 ) ,V ( 1 1 ) ,1 1 = 1 ,L PF ) 
CD TO P9R 

909 WPITF (NOT, 3001 ) 

WRITF ( NOT , 3002 ) LVlIl) 

WRITE (NOT, 3002) (LVC 1 1 ) , V ( 1 l ) ,1 1 = I ,LP2 E ) 
999 CALL ZZEOMP (5HYL0PD , NE PR DR ) 

END 
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Sl'fcFPlIT INE YM0D2A (NUTM ,NUTK ,NUTZ , W2 , W, FR E 0 , 

* NW ,\U,SHIFT,TCLZ ,T0LW2 ,MAXTT,I FPRNT, 

* V,LV»A,S,KVIN,KA*NUT1,NUT?,NUT3,NUT4 t NUT 5,NUT6) 
DIMFNSION V(l), LV(!), W2(!), W(l), FREQUI, A(KA,1), S(KA,1I 
DIMENSION IH(IO) 

DATA NI7 »NOT/S ,6/ 

CALCULATE MDDF SHAPES IPHJ) AND NATURAL FRFOUFNCIFS CF 
( —W?<- ‘ MAS S ) + (STIF ) )*(PHI I = (Cl USING 1TFRATIVF P AYLE IGH-R ITZ 
METHOD OF DR* JOHN ADMIRE. COMPOSITE STRUCTURE TECHNIQUE, 

NON— SWEEP INO VERSION, SPARSE PROGRAMMING LOGIC. 

THE MASS (NUTM I MATP IX SHOULD BE REAL, SYMMETRIC. 

THE STIF (NUTK) M ATR IX SHOULD BE REAL, SYMMETRIC. 

THE FIRST FLEMFNT OF FACH MODE SHAPE IS MADE POSITIVE. 

MODES APE NORMALIZED SUCH THAT (PHI IT* ( MASS 1 * (PH I ) = 1. 

CALLS FORMA SUBROUTINES .... 

BTAEA?,EIGN! A,INV4 ,MODEIX,NAME , PAGEHD , TIMCHK , KRI TE ,WRITIM, 
XLC'RD , YAABP ,vp SL3A, YDCM3A,YDTDS ,YIN ,YINI ,YLQRD ,YMULT1, 
YMUIT? ,YMULT^-» YNOZFR , YOUT ,YOim , YPART ,YRV1 ,YSTOD , YSYMLH, 
YSYMUH,VTRANr,YWRITE,ZZFCME. 

DEVELOPED FY RL KCHLFN, WA PENFIFLD, RA PHILIPPUS. MARCH 1972. 

LAST REVISION FY RL WOHLEN FOR NASA. MAY 1976. 

SUBROliT INt ARGUMENTS 

NUTM = INPUT LOGICAL NUMBER OF UTILI T Y TAPE CF MASS MATRIX. 

NUTK = INPUT LOGICAL NUMBER OF UTILITY TAPR CF STIF MATRIX. 

NUTZ = INPUT LOGICAL NUMBER OF UTILITY TAPF OF CALCULATED MODES. 

MAY BE USED TP TNPUT INITIAL RAYLEIGH VECTORS. 

W2 = OUTPUT VECTOR OF EIGENVALUES (OMEGA SQUARED). SIZE(NU). 

W = OUTPUT VECTTR OF CIRCULAR FR ECUFNCY (OMEGA). SIZF(NU). 

FREQ = OUTPUT VECTOR OF FREQUENCY (CMFGA/2P I I . SIZE ( NU ) . 

NW = NPUT NUMPFR OF MODES WANTED. ITERATIONS STOP WHEN NW 

CONS ECUT TVE MODES APOUND SHIFT POINT CONVERGE . 
HOWEVER, ALL NU MODES AND FREQUENCIES APE OUTPUT 
FOR LATER SELECTION. 

MU = INPUT NUMPFR OF MODES USED. MUST PE .GE. NW. 

SHIFT = INPUT SHIFT IN ( ST 1 F I -SHI FT ( M A SS ) . 

CONVERGENCE WILL PE ABOUT THIS VALUE. 

TOLZ = 1NPLO TOLERANCE ON ZERO W2. 

T0LW2 - INPUT CONVERGENCE TOLERANCE CN NON-ZERO W2. 

MAXIT = INPUT MAXIMUM NUMBER CF ITERATIONS. A GOOD VALUF IS 20. 
IEPRNT = INPUT = 1 HUNT INTFRME LI ATE P F SUL T S . 

V = INPUT VECTOR WORK SPACE. DIMENSION GREATER THAN * (N-l). 

LV - INPUT VECTOR WORK SPACE. DIMENSION GREATER THAN 6 * (N-l). 

A = INPUT MATRIX WORK SPACE. EQUIVALENCE TO LV AT KV/3-U. 

S = INPUT MATP IX WORK SPACE. EQUIVALENCE TO V AT KV/3+1. 

KVIN = I N r ‘ L'T DIMENSION SIZE CF V,LV IN CALLING PROGRAM. 

KA = INRUT ROW DIMENSION OF A,f IN CALLING PROGRAM. 

NUT1 = INPUT LOGICAL NUMEER OF UTILITY TAPE. 

NUT 2 = INPUT LOG I CAL NUMBER OF UTILITY TAPF. 

NUT? = INPUT logical numper of UTILITY T A P r • 

NUT4 - INPUT LOGICAL NUMBER OF UTILITY TAPE. 

NUTS = INPUT LOGICAL NUMPER OF UTILITY TAPE. 

NUTfr = INPUT LOGICAL NUMBER CF UTILITY TAPE. 
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C NFPROR EXPLANATION 

1 = NUMBER OF Mr*DES USED IS LESS THAN NUMBER WANTED. 

2 - NUMBER OF MODES USED FXCFFOS DIMENSION SIZE. 

3 = SIZE EXCEEDS KV. 


2010 FORMAT 

* 

* 

* 

* 

* 

* 

2020 FORMAT 
2030 FORMAT 

* 

* 

6000 FORMAT 
6020 FORMAT 


(/// 54X » 1 9H ( SUBROUTINE YMOD2A ) 

/// B2X » 19HN0. MODFS WANTED = 13, 

/// B2X » I9HND. MOOES USED = I?, 

/// 52X , 6HSHIFT = E15.8, 

/// 52X, 11FTOL ZERO = F15.E, 

/// 52X, l’HTOL W2 = EI5.B, 

/// E2X , BHMAXIT = 13) 

<// 7H ITFR = 13, 5H W2. / (1CX I0F1I.3) ) 

(/ 1?X 49HCCNVFRGENCE VALUES. EITHER W2CITER) IF (W2CITER) 
46H.LT. 2ER0 TOLER ANC F ) OR ( W2 f I TER )-W2 ( ITER-1 ))/ 
9HW2CITFR). / C10X 1CE11.3) ) 

( 2 OH NO. OF MODES BELOW ,E10.3,3H = ,IB) 

( 1 7X 16HITERAT10N TI?-1E = F1C.3,7H CP SEC F15.3,7H PP SEC) 


CALL TIMCHK (6HTBEGIN) 


KV = 4* C K VI N/4 ) 

CALL PAGFHD 

WRITE ) NOT, 2010 ) NW, NU, SHIFT ,TOLZ ,T0LW2 ,MAXIT 

IF (NU .LT. NW) GO TO 999 

IF CNU.CT.KA) CO TO 99Q 
P FWIND NU TM 

CALL YTNI (NUTM,IH, 1,10) 

N = IH(1) 

IF CKV.LT ,fc*(N-l ) ) GO TO 999 

IF (IFPRNT.FC.l) CALL YWPITF ( NUTM , 4HMA SS , V,L V,K V ) 
IF (JFRRNT.FO.I ) CALL YWPITF ( NUTK ,4HSTTF , V,L V,K V ) 


NE RRCR=1 
NEPP0P=2 


NE PR CR = 3 


CALCULATE IK-HAT) = (STIF) - SHIFT*CMASS ). 

CALL TIMCHK C 6HYAABB ) 

CALL YAAB ' ( 1 . , NITK ,-SHl F T , NUTM,NUT3 , V ,L V,K V,NUT4 t NUT5 ) 

CALL TIMCHK ( 6HYAABB ) 

IF (IFPPNT.FO.l) CALL YWPITF ( NUT3 , 5HK-HAT ,V,LV,KV) 

C DECOMPOSE (K-HAT) INTO ( U**T ) * ( -D-) * (U ) . 

CALL TIMCHK (6HYDCM3A) 

CALL YDCM3A (NUT3 ,NU7 1 ,NUT2, V ,LV,KV,NUT4, NUTS ) 

CALI TIMCHK ( 6HYDCM3A ) 

C CALCULATE NUMBER C’F NEGATIVE PLOTS ( FRFCUENC I FS) BELOW SHIFT POINT. 
REWIND NUT 2 

CALL YIN (NUT 2 ,V « 1 ,N) 

IF (IEPPNT .EC. 1) CALL WRITE { V, 1 , N , IHD, I ) 

NKCUNT = 0 
DO IE MKPUNT=1,N 

IF t VCMKOUNT ).LT.O. ) NKOUN T=NKCUNT + 1 
15 CONTINUE 

' WRITE (NPT,6C0C) SHIF T , NKOUN T 

C GENERATE INITIAL RAYLEIGH VECTORS. 


CALL TIMCHK ( 6HYPV1 ) 
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C 


C 

c 


c 


c 


CALL YR VI (NUTZ,N*NU,V t LV,KV,NUTA f NUT5,NUT6) Z 

CALL TIMCHK (6HYRV1 ) 

IF (IFPRNT.FC.il CALL YWRITF (NUTZ.6HZ-IN ,V,LV t KV) 

CALL TIMCHK (6HYMULTJ) 

CALL YMULT1 (NUT3,NUTZ,NUT4, V,LV,K V,NUT6 1 4 

CALL TIMCHK (6HYMULT1) 

IF (IFPRNT.EQ.l) CALL YWRITE ( NUT4, 2HKZ ,V, LV,K V) 

BEGIN ITERATION LOOP. 

CALL PAGEHD 
ITER = 0 

20 1TFR = ITEP+1 

CALL TIMCHK (6HYMULT21 

CALL YMULT? (NUTZ ,NUTA t NUT3» A * S «V « L V,KV»K A ,NUT6) 3 

CALL TIMCHK (6HYMULT2 ) 

CALL TIMCHK (6HYMULT1 ) 

CALL YMULT! <NUTm t NUTZ *NUT5» V,LV»KV.NUT6 ) 5 

CALL TIMCHK ( 6HYMULT1 ) 

IF (IFPRNT,EC.l ) CALL YWRITE ( NUT5 , 2HMZ t V, LV,K V) 

CALL TIMCHK (6HYMULT2) 

CALL YMUL 7 2 (NUTZ ,NUT 5 *NUT4, A, S, V t L V,KV ,KA *NUTfc) 4 

CALL TIMCHK (6HYMULT2 ) 

CALL TIMCHK (6HYST00 ) 

CALL YSTPP (NUT4 t A ,NR , NC,KA ,K A, V , L V,KV ,NUT6 ) A 

CALL YE T OD INUT3 »S»NR ,NC »KA »KA»V»LV»KV*NUT6) S 

CALL TIMCHK (6HYST0D » 

IF (IFPPNT .FC. I) CALL WRITE ( A,NU ,NU, AHMRAP ,KA ) 

IF (IFPPNT .FQ. I) CALL WRITE (S ,NU .NU» AHKRA* ,KA 1 

CALL TIMCHK (6HM0DE1X) 

CALL MODE IX ( A ,S * W2 ,NU ,T0LW2 ,K A ) A 

CALL TIMCHK ( 6HM0DE1 X ) 

IF IIFPRNT .EQ. 1) CALL WRITE (A, NU,NU, 2HY * , KA) 

CALL TIMCHK (6HCVTEST ) 

UNSHIFT W2. 

DP 25 J=1 »NU 

V( J) = AE S(W2( J) ) 

LV(J) = J 

25 W2(J> = W2(J) ♦ SHIFT 

WRITE (NET, 2020) ITF R , ( W2 ( J ) , J = 1 ,NU ) 

IF ( TTF R .EC. 1 ) GP Tp feO 

STORE CONVERGENCE VALUES OF W 2 IN W. < LAST ITER CF W2 I- IN FREQ t • 

DC 28 JsIfNU 
W( J) = W2 ( J) 

IF ( ARE (W2 ( J ) 1 .GT. TOLZ) W(J1 = ( W2 C JJ-FR FQ I J 1 >/W2 ( J > 

28 CPNTINUF 

WR ITF (NCT * 2030) (W ( J ) , J= 1 ,NU ) 

IF (ITER.CE .MAXIT) GO TP 7C 

FIND START AND FND LCCATIPN OF f W2-SHI FT ) OF BAND WIDTH NW ABOUT SHIFT 
DP 37 J=I,NW 
IMIN = J 
VMTN ~ V(J) 

LMIN = IVU) 

DP 36 I=J,NU 

IF (VMIN .IF. V( I )) GO TP 36 
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IMIN = I 
VMIN = VCI) 

LMIN = LVCI) 

36 CONTINUE 
V(IMIN) = V(J) 

LV(IMIN) = LV(J) 

VCJ) = VM IN' 

37 LVCJ) = LMIN 
JS = LVCI) 

je = lvci > 

DC 38 J = 1 ,NW 

IF CLVCJ) .LT. JS) JS=LV( J) 

IF CLVCJ) .GT. JF) JE=LV(J) 

38 CONTINUE 

C TEST W2 FCF CONVEFGE NCE OF NW CONSECUTIVE MODES APOUT SHIFT POINT. 
DO 45 J = JS » JE 

IF ( ABSCW2C J) ) .LT. TCLZ) CO TO 45 
IF CAP.SCWCJ)) .GT. TOLW2 ) GO TO 50 
45 CONTINUE 
GO TO 70 
50 CONTINUE 

C STORE LAST ITERATION VALUE OF W2 IN FRFQ. 

60 DO 62 J = 1 ,NU 
62 FREQCJ) = W2CJ) 

CALL TIMCHK 

C IMPROVE RAYLEIGH VECTOPS. 

CALL TIMCHK 

CALL YMUL T4 (NUT 5 ,A ,NUT4 ,S ,V , LV,KV, KA ,NUT6 ) 

CALL TIMCHK 

IF CIFPRNT.EO.l ) CALL YWRITE (NUT4 , IHG, V, L V,K V) 

CALL TIMCHK 

CALL YPSL3A (NUT 1 ,NUT 2 ,NUT4, NUTZ , V ,LV, KV , NUT 5 ,Nl»T6 ) 

CALL TIMCHK 

IF CIFPRNT.FO.I ) CALL YWRITE ( NUT2 , 1HZ, V, LV,K V ) 

GO TO 20 

C END ITERATION LOOP. 

C 

C GET W,FREO,MODE SHAPES. MAKE FIRST ELEMENT OF EACH MODE 
C SAVE ALL MODFS AND FREQUENCIES USED CNU) FOR LATER SELECTION. 

70 DO 72 1=1 ,NU 

WCI) = SCR T CAESCW2CI))) 

72 FREO(I)= .15915494 * WCI) 


C ALL 

CALL YMULT4 (NUT Z ,A ,NUT 1 ,S ,V , L V,K V, KA ,NUT6 ) 

TIMCHK 

C6HYMULT4) 

1 

CALL 

TIMCHK 

(6HYMULT4) 


CALL 

TIMCHK 

(6HMD1P0S) 



RFWIND NUT1 
REWIND NUT 7. 

IVSMI = KV-N 
IVS = IVSMI ♦! 

DO 73 1= TVS ,KV 
VCI) = 1. 

CALL Y IN 1 (NUT 1 ,IH, 1 ,10) 
CALL YOUTI (NUTZ ,IH, 1,10) 
NGZ = IHC 3 ) 


(6HCVTEST) 

C6HYMULT4) 

u 

C6HYMULT4) 

(6HYBSL3A) 

Z 

(6HYBSL3A) 


POSITIVE. 


73 
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DO 75 IGZ=l t NGZ 

CALL YI> ‘I tNUTl »IH» 1*10 

NELGZ = IHC1) 

CALL YINI (NUTl,LV,l, NELGZ) 

CALL YIN tNl'TI ,V*1* NELGZ) 

DO 74 1=1, NELGZ 
IZ = LV ( I )/!OPCOO 
JZ = LV(I)-10000ff*IZ 

IF (IZ.FC.l .AND. V( 1 ) .LT. 0. ) V< IVSM1+JZ)=-1. 

74 Vfl) = VII )*V( IVSM1+J Z ) 

CALL YOUTI CNUTZ,IH,1,1C) 

CALL YOUTI !NUTZ,LV,1, NELGZ) 

75 CALL YOUT (NL’TZ ,V*1, NELGZ) 

CALL TIMCHK I6HMD1P0S) 
CALL TIMCHK (6HTPRINT) 

RETURN 

999 CALL ZZBOMB (6HYMCDE2 ,NERRCR ) 

END 
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SUBROUTINE YMULT (NUTA *NUTB , NUTZ ,V*LV»KV,NUT1 I 
DIMENSION V(1),LV(1), MHEAD ( 10 ) 

DATA NITtNCT/5,6/ 

C 

C SPARSE MATRIX MULTIPLICATION (A * B = Z). 

C KV/4 MUST BF ECUAL TO OR GREATER THAN, 

C Cl I NUMBER OF COLUMNS OF MATRIX A 

C AND < 2 ) NUMBER OF COLUMNS OF MATRIX B. 

C CALLS FORMA SUBROUTINES YIN ,YINI » YLORD ,YNOZER,YOUT ,YOUTI , 

C YPART ,YSYMLH,YSYMUH,ZZBOMB* 

C DEVELOPED BY R A PHILIPPUS. APRIL 1969. 

C LAST REVISION BY RL VfclHLEN FOR NASA. MAY 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

NUTB = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX E IS STORED. 

NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

NERROR EXPLANATION 

1 = SIZE LIMITATION EXCEEDED. 

2 = INCOMPATIBLE MATRICES. 

NPARTZ=0 
NNZZ=0 
NRFC=0 

CALL YLORD (NUT A ,V, L V ,K V,NUT1 ,NUTZ ) 

CALL YLORD (NUTB ,V, LV ,KV, NUT1 , NUTZ ) 

C GET (A) HEADER INFORMATION. 

REWIND NUTA 

CALL YIN1 (NUTA »MHE AD, 1,10) 

NR A = MHFAD(l) 

NCA = MHF AD ( 2 ) 

NNZ A = MHEAD (A ) 

ISHAP = MHF AD ( 7 ) 

NPCT A= 1C0TNNZA/NR A/NC A 
IXYZ3=ISHAP 

IF (ISHAP. FG.5HWHOLB -OR. ISHAP. E0.4HDI AG ) GO TO b 

IF ( ISHAP .EC .SHLOWER ) CALL YSYMUH ( NUTA, V , LV, KV, NUT1 , NUTZ ) 

IF (ISHAP. F0.5HUPPER) CALL YSYMLH CNUTA,V , LV, KV,NUT1 ,NUTZ ) 

REWIND NUTA 

CALL YINI (NUT A »MHE AD ,1 , 10 ) 

NNZ A = MH EAD (4 ) 

NPCT A=1 C0*NNZA/NP A/NC A 
b NPARTA = MHEAD (3 ) 

C GET (B) HEADER INFORMATION. 

REWIND NUTB 

CALL YINI (NUTB ,MHEAD, 1,10) 

NRB = MHEAD (1 ) 

NCB = MHEAD (2) 

NNZP = MHEAD (4) 

ISHAP = MHE AD( 7) 
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C 


C 


10 


NPCTB=1P0*NNZB/NR6/NCB 
IXYZ4= ISHAP 

IF (ISKAP.E0.5HWHCLE .OR. ISHAP.EQ.4HDIAG > GO TO 10 
IF CISHAP.E0.5MLOWER) CALL YSYMUH (NUTB ,V * LV» KV»NUT1 ,NUTZ ) 

IF CISHAP.EQ.5HUPPER) CALL YSYMLH CNUTB ,V, LV,KV,NUT1 ,NUTZ) 
REWIND NUTB 

CALL YINI CNUTB, MHE AD* I *10) 

NNZB = MHEADC4) 

NPCTB=1CC*NNZB/NRB/NCB 

IF (NNZA.FO.O .OP. NNZB.EQ.OI GO TO 70 

NPARTB = MHEADC3 ) 

NERROR*! 


IF (NCA.GT.KV/4 .CR. NCB.GT.KV/4) GO TO 999 
IF (f RB.NE.NCA ) GO TO 999 


N ERRORS 


IZ=0 

LPBS=KV/4+l 
LPEE=LPBS-1 
LCS=KV/2+l 
LCE=LC$— 1 +NC6 
LCC$=LCE+1 
LCCF=LCE 
NNZ=KV— LCCS+1 
REWIND NUT1 

DO 15 I=LCS ,LCE 
15 V(I)=0. 


LOOP ON ( A 9 PARTIONS. 

DO 55 1=1 .NPARTA 
C GET (A) PARTITION INFORMATION. 

CALL YINI (NUT A ,MHE AD, 1 * 10 ) 

NNZPA = MHEAD(l) 

LFFLPA = MHEADC2) 

LLELPA = MHEA0C3) 

CALL YINI (NUTA,LV, 1, NNZPA) 

CALL YIN <NUTA,V,1 , NNZPA) 

K=LPPS 
ITRPL=0 
REWIND NUTB 

CALL YINI CNUTB ,MHE AD, 10, 10) 

NREAD=0 

C 

DO 50 INA=1, NNZPA 
I A=LV C INA )/l OOOCO 
JA=LV( INA )-l 00000*1 A 

IF (IA.FQ.IZ .AND. ITRBL.EQ. 1 ) GO TO 50 
ITRBL=0 

IF CIA.FO.IZ) GO TO 30 
REWIND NUTB 

CALL YINI (NUTB ,MHE AD, 10, 10) 

NREAD=0 

C 


DO 25 INC = LCS,LCE 
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IF ( V« IMG I.EQ.O.) GO TO 25 
LCCE=LCCE I 

V ( LCCE ) =V ( INC) 

LVf LCCE >=I2Z+lNC-KV/2 
V( INC )=0. 

IF <LCCE.LT.KV) GO TO 25 
CALL YOUTI (NUT1 ,LV, LCCS,LCCE ) 

CALL YOUT (NUT1,V*LCCS»LCCE) 
NREC=NREC+1 
NNZZ*NNZZ-»NNZ 
LCCE=LCE 
25 CONTINUE 

1 2= I A 

IZZslOOCOO*IZ 

K=LPP« 

30 IF CK.LE.LFBE .AND. NREAD.GT.O) GO TO 40 
35 IF (NRF AD.EO.NPARTB) ITRBL=1 
IF IITRBL.EQ.1J GO TC' 50 
CALL YINI (NUTB *MHE AD»1 »!0 ) 

NNZPP s MHFADd ) 

LFELPB = MHEAD <2 ) 

LLELPB = MHEAD13) 

LPBE=LPBS-1+NNZPB 

CALL YINI <NUTB,LV,LPBS*LPBE) 

CALL YIN <NUT6,V,LPBS f LPBE) 

NREAD=NREAD+1 

K=LPB? 

40 DO 45 INB=K tLPBE 
K = INB 

IB=LV<INB)/1 00000 
IF (Ifi.GT.JA) GO TO 5 0 
IF flB.LT.JA) GO 70 45 
JBZ*LV< INB)-10000C*IB 
INZ=KV/ *:♦ JBZ 

V I INZ ) »V ( INZ ) +V{ INA ) *V ( INB ) 

45 CONTINUE 

GO TO 35 
50 CONTINUF 
55 CONTINUE 

DO 60 1=LCS,LCF 

IF IVIII.FO.O.) GO TO 60 

LCCE=LCCE+I 

V < LCCE I - V < 1 1 
LV(LCCF) = IZZ + I~KV/2 

IF (LCCE.LT.KV) GO TO 60 
CALL YOUTI (NUT 1 *1. V t LCC S * LCCE ) 

CALL YOUT <NUT1,V,LCCS»LCC£) 
NRFC=NREC+1 
LCCE=LCE 
NNZZ=NN2Z+NNZ 
60 CONTINUE 
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IF (LCCE.EQ.LCE) GO TO 7C 
NNZ=LCCE— LCCS+1 
NNZZ=NNZZ +NNZ 
NREC=NREC+1 

CALL YOUTI (NUT1 »LV* LCCS *LCCE ) 
CALL YOUT (NUT1,V,LCCS,LCCE) 
70 REWIND NUTZ 
MHEAD(l) = NRA 
MHEADJ2 J = NCR 
MHEADt 3 ) = NREC 
MHEAD(4) = NNZZ 
MHEAD(5) = 5H0PDER 
MHEADC6) = 0 
MHE AD ( 7 ) = 5HWH0LE 
CALL YOUTI (NUTZ *MHE AD, 1,10) 

IF (NNZZ.GT.O) GO TO 75 
DO 72 1=1,10 
72 MHEAD(I) = 0 

CALL YOUTI (NUTZ ,MHE AD, 1 , 10 ) 
CALL YOUTI CNUTZ ,MHEAD*1*2 ) 
CALL YOUT (NUTZ, V,I ,2) 
RETURN 
75 LZE=KV 

REWIND NUT1 

DO 100 1=1, NREC 

IF (I.EQ.NPEC) LZE=LCCS— I+NNZ 

NNZP = LZE— LCCS+1 

CALL YINI (NUT! ,LV,LCCS*LZE ) 

CALL YIN (NUT1 »V»LCCS»LZE ) 

MHEAD ( 1 ) = NNZP 

MHEAD ( 2 ) = LV(LCCS) 

MHEAD( 3 ) = LV(LZE) 

CALL YOUTI (NUTZ, MHE AD, 1,10) 
CALL YOUT: (NUTZ, LV, LCCS, LZE ) 

100 CALL YOlfT (NUTZ ,V, LCCS , LZE ) 

CALL Y PART f ,'IUTZ ,V,L V ,KV, NUTl ) 
RETURN 

999 CALL 2.ZBCME ( 5HYMULT ,NERROR ) 
END 
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SUBROUTINE YMULTI (NUTA, NUTB , NUTZ, V,LV,KV,NUTI ) 

DIMENSION V(1),LV(1),MHEAD(10) 

DATA NIT, NOT/5,6/ 

C 

C SPECIAL SPARSE MATRIX MULTIPLICATION I A * B = Z). 

C B AND Z ARE DENSE MATRICES. 

C KV/A MUST BE EQUAL TO OR GREATER THAN, 

C (1) NUMBER OE COLUMNS OF MATRIX A 

C AND (2) NUMBER OF COLUMNS CE MATRIX B. 

C CALLS FORMA SUBROUTINES YIN ,YINI , YLORD ,YNOZER,YOUT ,YOUTI , 

C YPART ,YSYMLH,YSYMUH,ZZBOMB. 

C DEVELOPED BY R A PHILIPPUS. AUGUST 1972. 

C LAST REVISION BY WA BENE I ELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS ( ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

C NUTB = LOGICAL NUMBER OE UTILITY TAPE ON WHICH MATRIX R IS STORED. 

C NUTZ = LOGICAL NUMBER OE UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NEPROR EXPLANATION 

C 1 = SIZE LIMITATION EXCEEDED. 

C 2s INCOMPATIBLE MATRICES. 

KVCA = KV/A 

CALL YLORD (NUT A ,V , L V ,KV ,NUT1 ,NUT2 ) 

CALL YLORD (NUTB ,V, LV,KV,NUT1 , NUTZ ) 

REWIND NL’TA 

CALL YINI (NUT A ,MHE AD,1 , 10 ) 

NR A » MHF AD ( 1 ) 

NCA = MHEAD( 2 ) 

NNZ A = MHEAD ( A) 

ISHAP = MHEAD( 7) 

IF ( ISHAP. EQ.5HWHCLE .OR. ISHAP. FQ. AHDIAG ) GO TO 5 

IE ( ISHAP. EC. 5HL0WER) CALL YSYMUH (NUTA, V, LV,KV,NUT1 , NUTZ ) 

IF ( ISHAP -E0.5HUPPER) CALL YSYMLH ( NUTA ,V , LV, KV»NUT1 ,WU TZ) 

REWIND NUTA 

CALL YINI (NUTA »MHE AD, I ,10 ) 

NNZ A = MHEAD(A) 

5 NPARTA = MH E AD ( 3 ) 

REWIND NUTE 

CALL YINI (NUTR,MHEAD,1 ,10) 

NRB = MHF AD ( 1 ) 

NCR = MHE AD (2 ) 

NPARTB = MHE AD (3 ) 

NNZB = M.HEAD(A) 

NE RR0R = 1 

IF (NCA.GT.KV/A NCB.GT.KV/A) GO TO 999 

NERR0R=2 

IE (NRB.NE.NCA) GO TO 999 


C 


1Z=C 
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LBS=KV04+1 
LZS = 2*KV04*1 
NNZZ = NNZB 
NPARTZ = NPARTB 
NNZPZ = KV04/NCB*NCB 
MHFAD(l) = NR A 
MHEAD (2 ) = NCB 
MHFADC?) = NPARTZ 
MHEAD(4) = NNZZ 
MHEAD( 5 I = 5HCP.DE R. 

MHEAD(6 ) = 4*KV04 
MHEAD(7) = 5HWH0LE 
REWIND NUTZ 

CALL YOUTI (NUTZ,MHPAD,1,10) 

READ A ONE TIME , EACH PARTITION AS REOU1RED. 

CALL YINI (NUTA.MHEAD, 1,10 

NNZPA = MHEAD(l) 

LFELPA = MHEAD ( 2 ) 

LLELPA = MHEAD <3 ) 

NNZARD - NNZPA 

CALL YINI (NUTAtLV, 1 ,NNZPA) 

CALL YIN (NUT A * V,l, NNZPA) 

INA = 0 
NFRPZ = 1 

DC 55 IPARTZ=1 , NPARTZ 

IF (IPARTZ.FQ. NPARTZ) NNZPZ=NNZZ-( I PARTZ-1 )*NNZPZ 
NRPZ = NNZPZ /NCB 
NLRP?. = NFRPZ-l+NRPZ 
L2E -- LZS-1+NN2PZ 
INZ = LZS— 1 
DO 20 I=NFRPZ,NLRPZ 
DC 2C J-ltNCe 
INZ = INZ+I 
V(INZ) = 0. 

20 LV(INZ) = 100000*I+J 

21 REWIND N'JTB 

CALL YINI (NUTB,MHEAD f l ,10) 

C 

C READ ALL CF P FOR EACH PARTITION OF A OP Z. 

DO 50 IPAPTB = ! ,NPARTB 

CALL YINI (NUTB, MHEAD *1,10) 

NNZPB = MHEAD ( 1 ) 

LFELPB = MHEAD ( 2 ) 

LLELPB - MHEAD (3 ) 

LBE = LPS-1+NN2PB 

CALL YINI (NUTB *LV, LBS ,L£E ) 

CALL YIN (NUTB, V, LBS, LBE) 

INA = 0 

NFRPB = LFPLPB/1 CCCOC 
NLRPE = LLELPB/1 00000 


25 IF ( INA. LT. NNZPA) GO TO 30 

IF (IPARTB.LT.NPARTB) GO TO 50 
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IF (NNZARD.EQ.NNZA) GO TO 50 
CALL YIN I INUTA »MHE AD» 1 » 1C ) 

NNZPA = MHEADU ) 

LFELPA = MHEAO (2 ) 

LLELPA = MHEAD (3 ) 

CALL YINI (NUTA,LV,1, NNZPA) 

CALL YIN (NUTA, V,l, NNZPA) 

NNZARD = NNZARD+NNZPA 
INA = 0 
GO TO 21 
C 

30 ’NA = INA +1 

IA = LV ( INA ) /l OOOCO 
IF (IA.LE.NLRPZ) GO TO 35 
GO TO 50 
C 

35 IF (IA.LT.NFRPZ) GO TO 25 
JA = LV ( INA ) -1C0000*I A 
IF (JA.LT.NFRPB) GO TO 25 
IF (JA.GT.NLRPB) GO TO 25 
LZ = ( I A— NFPPZ ) *NCB+L Z S— 1 
LB = < JA-NFRPB)*NCB+LBS-1 
IF (LB+NCB.GT.LBE ) GO TO 50 
C 

DO 40 JZ=1,NCB 
LZ = LZ+1 
LB = LP+I 

40 V ( LZ ) = V (LZ)+V( INA)*V(LB) 

GO TO 25 

C 

50 CONTINUE 

NFRPZ = NLRPZ+1 
MHEAO (2) = NNZPZ 
MHEADf 2 ) = LV(LZS ) 

MHFADC3) = LVCLZE) 

CALL YOUTI (NUTZ »MHE AD, 1 1 10 ) 
CALL YOUTI (NUTZ ,LV,LZS,LZE ) 

55 CALL YOUT (NUTZ, V,LZS,LZE) 
CALL YPART (NUT2 ,V, L V,KV,NUT1 ) 
RETURN 
C 

999 CALL Z2B0MB (6NYMULT1 ,NE RROR ) 
END 
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SUBROUTINE YMULT2 (NUT A, NUTP ,NU7Z , A ,e ,V ,L V t KV ,KR « NUT1) 

DIMENSION V(l), LV(1), A(KR,1), B(KR,1), MHEAD(IO) 

COMMON / LWRKV1 / W(500) 

DATA NIT*NCT/5,6/ 

C 

C SPECIAL SPARSE MATRIX MULTIPLICATION ( A**T )* ( B) = ( 2) . 

C A»B,Z ARE DENSE MATRICES. 

C 2 IS SYMMETRIC. 

C KV/A MUST BE FOUAL TO OR GREATER THAN, 

C (11 NUMBER OF COLUMNS OF MATRIX A 

C AND (2) NUMBER OF COLUMNS OF MATRIX B. 

C CALLS FORMA SUBROUTINES YDTOS ,YIN ,YINI , YLORD ,YNOZER,YGUT , 

C YOUTI , YPAPT ,YSYMLH , YSYMUH>ZZBOMB. 

C DEVELOPED BY R A PHILIPPUS. AUGUST 1972. 

C LAST REVISION BY WA BENF IELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA = LOGICAL NUMPE* OF UTILITY TAPE ON WHICH MATRIX A IS STORED- 

C NUTB = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX B IS STORED. 

C NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

C A = MATRIX WORK SPACE. 

C B = MATRIX WORK SPACE. 

C V = VFCTOP WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C KR = ROW DIMENSION OF A,B IN CALLING PROGRAM. 

: NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C NEPROR EXPLANATION 

C 1 = SI Z F LIMITATION EXCEEDED. 

C 2 - INCOMPATIBLE MATRICES. 

C 3 = INCOMPATIBLE MATRICES. 

C A - INCOMPATIBLE MATRICES. 

C 5 = INCOMPATIBLE MATRICES. 

C 6 = MORE THAN 500 COLUMNS IN MATRIX B. 

C 

CALL YPART (NUTA ,V « L V»K V,NUT1 ) 

REWIND NUTA 

CALL YINI (NUTA ,MHE AD, 1 , 10 ) 

NR A = MHEAD(l) 

NCA = MHFAD ( 2 ) 

NNZA * MHEAD(A) 

NP ART A = MHEAD ( 3 ) 

REWIND NUTP 

CALL YINI (NUTB ,MHE AD, 1 , 10 ) 

NRB = MHEAD ( 1 ) 

NCB = MHEAD ( 2 1 
NNZB = MHEAD (A) 

NPARTB = MHEAD ( 3 ) 

NERRCRsl 

IE (NCA.GT.KV/A .OR. NCB.GT.KV/4J GO TO 999 

NERROR=2 

IF (NRB.NF.NPA5 C-C 70 999 
C 


NE RR0R=3 
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IF (NPARYA.NE »NP AR1 b ) GO TO 999 

IF (NNZA.NE.,,NZB ) GO TO 999 

LZS = KV/2+1 

L2E = LZS-1+NCA*NCA 

REWIND NUT1 

DC 30 IPART=1 tNPARTA 

CALL YINI (NUTA,MHEAD,i,10) 

NNZPA = KHEAD(I) 

IRF = MHEAD ( 2 ) /I 00000 
IRl - MHEAD(3)/10000C 
NR * IRL-IRF+1 
CALL VIN (NUT A ,V, 1 , 1 ) 

CALL YIN (NUTA,V,1 , NNZPA) 

f ALL YINI (NUTB ,MHE AD, 1 , 10 ) 

IF (MHE AD (l).NE. NNZPA) GO fO 999 

LPBS = NNZPA+1 

LPBE = NNZPA+MHE AD( 1 ) 

CALL YIN (NUTB ,V,L PBS , LPP S ) 

CALL YIN (NUTB»V, LPBS, LPBE) 


NERF OR -4 


NFRRCR=5 


STATFMENTS FROM ATXBB2 

NERR0R=6 

IF (NCA.GT.5C0) GO TO 999 
DO 20 J=1 ,NCA 
DO 10 1=1 ,J 
W(I) = 0. 

DO 10 K=1 ,NR 
IN A = (K-l )*NCA*1 
INB = NNZPA-MK-1 )*NCA+J 
10 W(I) = W( I ’'♦VI INA )*V ( INB ) 

DO 20 1 = 1, J 
A ( I , J) = W ( I 1 
20 A ( J , I ) = W(I) 

30 CALL YOUT (NUT1,A,1, KR*NCA ) 

IF (NPAPT A.LF. 1 ) GO TO 40 

REWIND NUT1 

DO 35 IPART=2,NPAPTA 

CALL YIN (NUT1 ,B,1 , KR*NCA> 

DO 35 1=1, NCA 
DO 35 J=1 ,NCA 
35 A ( I , J ) = A ( I *J ) + B (I , J ) 

40 CALL YD TO S ( A ,NUTZ,NC A,NC A,KR *KR , V,LV,KV ,NUT 1 ) 

RETURN 

999 TALL ZZBOMB (6HYMULT2 ,NERROR ) 

END 
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SUBROUTINE YMULT4 'NUTA, B, NUTZ ,A , V , LV ,KV, KP ,NUTI ) 

DIMENSION VCI)» LV(1), A(KR,1), e(KR,l), MHFAD <10 1 
COMMON / LWRKV1 / W<500> 

DATA NlT, NOT/5, 6/ 

C 

C SPECIAL SPARSE MATRIX MULTIPLICATION <A * B = Z). 

C A,B,Z ARE DENSE MATRICES. 

C KV/4 MUST PE EQUAL TO OR GREATER THAN NUMBER OF COLUMNS CF MATRIX A. 
C CALLS FORMA SUBRf JTINES YIN ,YINI ,YLORO ,YOUT ,YOUTI ,YPART , 

C ZZBOMB. 

C DEVELOPED BY R A PHILIPPUS. AUGUST 1972. 

C LAST REVISION BY WA BENF I ELD . MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS <ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

C B = MATRIX. SIZE <NC A,NCA). 

C NUTZ = LOGICAL NUMEER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

C A = MATRIX WORK SPACE. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C KR = ROW DIMENSION OF A,E IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NFPROR EXPLANATION 

C 1 = SIZE LIMITATION EXCEEDED. 


REWIND NUTA 

CALL YT VI (NUT, MHE AD, 1 , 1G ) 

NR A = MHF AD < 1 ) 

NC A = MHE AD ( 2) 

NNZ A = MHE AD (4 ) 

NPARTA = MHEAD ( 3 ) 

NRB = NCA 
NCP = NCA 

IE (NCA.GT.KV/4) GO TO 999 
C 

REWIND NUTZ 

CALL YOUT! (NUTZ ,MHE AD,1 , IC ) 
DO 50 IPART=1, NPARTA 
CALL YINI (NUT A ,MHE AD,1 , 1C ) 

NNZ P A = MHEAD(l) 

IRE = MHEAD(?)/1CCCC0 
IRL = MHEAD(3)/I 00000 
NR = IRL-IRF+l 

CALL YTNT (NUT A ,LV, 1 ,NNZPA ) 

CALL YIN (NUTA,V,I ,NNZPA) 

DC AC 1=1 , NR 
DO 20 K=?,NCA 
INA = ( I— 1 ) *NC A+K 

20 W(K) = V ( INA) 

DO 40 J=I,NCA 
5*0. 

DO 30 K=1,NCA 


NEPR0R=1 
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30 S = S+W(KI*P(K,J» 

IN? & (I-D+NCA + J 
40 VC1NZ) = ? 

CALL YOUTI (NUTZ .WHEAO, 1,101 
CAr YOUTI (NUTZ »LV »1»NNZPAI 
50 CAu*. YOUT (NUTZ t V«1«NN2PA) 
RETURN 

999 CALL ZZBOMB (6HYMULT4»NERRCR ) 
END 
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SUBROUTINE YNOZER (NUTAZ t V,LV,KV,NUT1 ) 

DIMENSION V(1),LV(1), MHEADf 10) 

DATA NIT tNOT/5 *6/ 

REMOVE ZERO ELEMENTS FROM SPARSE MATRIX A TO GET SPARSE MATRIX Z. 
MATRIX A IS REPLACED BY MATRIX Z ON NUTAZ. 

CALLS FORMA SUBROUTINES YIN ,YINI ,YCUT ,YGUTI ,YPART ,ZZBOMB. 
DEVELOPED BY R A PHILIPPUS. OCTOBER 1968. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTAZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRICES A,Z ARE 
STORED. MATRIX A IS REPLACED BY MATRIX Z. 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

NERROR EXPLANATION 
1 = SIZE LIMITATION EXCEEDED. 

REWIND NUTAZ 
REWIND NUT1 
NREC=0 

CALL YINI (NUTAZ,MHEAD,1,1C) 

NRA = MHFAD(l) 

NCA = MHEAO ( ? ) 

NPART = VHE AD ( 3 ) 

NNZ A = MHFAD (4 ) 

MCKORD = MHE AD (5 ) 

MASHAP = MHEAD (7 ) 

IF (NNZA.EQ.O) GO TO 25 

NNZCK-NNZA 

NNZA=0 

DO 10 I =1 » NPART 

CALL YINI (NUTAZ*MHEAD,1,I0) 

NNZP = MHEAD ( I ) 

NERR0R=1 

IF (NN2P.GT.KV) GO TO 999 

IF (NN2P.GT .0 ) GO TO 5 

CALL YINI (NUTAZ, MHEAD, 1C, 10) 

CALL YINI (NUTAZ, MHEAD*10, 10) 

MHFAD(IC) = 0 
GO TO JO 

5 CALL YINI (NUTAZ, LV,1, NNZP) 

CALL YIN (NUTAZ, V*l, NNZP) 

NNZZ=0 

DO 8 J= 1 »NN2P 

IF (V(J I.EC.C.) GO TO 8 

NNZZ=NNZZ+1 

V ( MN 2 Z ) =V ( J ) 

LV(NNZZ )=LV( J) 

8 CONTINUE 
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£ 

IF (NNZZ.EO.C) GP TO 10 
NNZA=NNZA+NNZZ 
NRFC=NRFC+1 
MHEADID = NNZ2 
MHFADt 2 ) = LVII ) 

MHFAD ( 3 ) = LV(NNZZ) 

CALL YOL’TT (NUTI ,MHEAD, 1 ,10 J 
CALL YOUTi (NUTI ,LV,1 ,NNZZ) 

CALL YOUT (NUTI ,V,1,NNZZ) 

10 CONTINUE 

C 

IF (NNZA.EQ.NNZCK) GO TO 25 

REWIND NLITAZ 

REWIND NUTI 

MHEADtl) = NRA 

MHFAD ( 2 ) = NCA 

MHEAD(3) - NREC 

MHFAD (4) = NNZ A 

MHEADC5) = MCKORD 

MHEADl 7 > = MASHAP 

CALL YOUT I (NUTAZ, MHFAD, 1,10 

IF (NRFC.GT.O) GO TO 15 

MHEADtl) = 0 

MHFAD (2) = 0 

MHFAD 1 3 ) = C 

CALL YOUT I CNUTAZ, MHFAD, 1 ,10 J 
CALL YOUT I (NUT A 2, MHFAD, 1,2 » 

CALL YOUTT (Nl>TAZ,MHEAD, 1 ,2 ) 
RETURN 
C 

15 DO 2 C 1=1, NREC 

CALL YINI (NUTI , MHFAD, 1,10 

CALL YINI (NUTI ,LV, 1 , MHEAO ( 11) 
CALL YIN (NUTI, V,1 , MHEADtl )1 
CALL YOUTI (NUTAZ, MHFAD, 1,10 
CALL YOUTI (MJTAZ, LV,1, MHEADtl )) 
20 CALL YOUT (NUTAZ,V , 1 , MHFAD (11) 

C 

25 CALL YPART (NUTAZ, V, LV,KV, NUTI) 
RFTURN 
C 

999 CALI ZZBCKB (6HYN0ZER ,NERRCR ) 

END 
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SUBROUTINE YOUT (NUT,A,NS,NE ) 

DIMENSION AC1) 

WRITE OUT BINARY DATA FROM CORE AREA A ONTO PERIPHERAL UNIT NUT. 
CALLS FORMA SUBROUTINE ZZBCMB. 

DEVELOPED BY R A PHILIPPUS. NOVEMBER 1971. 

LAST REVISION BY WA BENFIELD FOR NASA. MAY 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUT = LOGICAL NUMBER OF UTILITY TAPE. 

A = VECTOR TO BE WRITTEN. 

NS = START LOCATION IN VECTOR A. 

NE = END LOCATION IN VECTOR A. 

NfRROR EXPLANATION 

1 = START LOCATION GREATER THAN END LOCATION. 

2 = END OF FILE ENCOUNTERED. 

3 = WPITF PARITY ERROR. 

NERRORM 

IF (NS .LE. 0 .OR. NS .G1 . NE) RETURN 

NERROR=2 

WRITE (NUT* ERR=998»END=999) (A ( I ) * I=NS*NE ) 

RETURN 
C 

998 

999 CALL ZZBOMB (4HY0UT *NERRC'R ) 

END 


NERRCR=3 
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SUBROUTINE YOUTI (NUT , IA ,NS ,NE ) 

DIMENS ION T A ( 1 ) 

C WRITE OUT BINARY DATA FROM CORF AREA IA ONTO PERIPHERAL UNIT NUT. 
C CALLS FORMA SUBROUTINE ZZBOM6. 

C DEVELOPED BY R A PHILIPPUS. MAY 1973. 

C LAST REVISION EY WA EENFIELD FOR NASA. MAY 1976. 

C 

C SUBROUTINE ARCUMENTS (ALL INPUT) 

C NUT = LOGICAL NUMBER OF UTILITY TAPE. 

C IA = VFCTOR TO BE WRITTEN. 

C NS = START LCCATICN IN VECTOR IA. 

C NE = END LOCATION IN VECTOR IA. 

C 

C NERROR EXPLANATION 

C 1 = START LOCATION GREATER THAN END LOC ’TON. 

C 2 = END OF FILE ENCOUNTERED. 

C 3 = WRITE PARITY ERRCR. 

C 

NERROR=l 

IF (NS .LB. 0 .OR. NS .GT. NE) RETURN 

NERRCR=2 

WRITE (NUT,ERR=998,END=999) ( I A ( I ), I=NS ,NE ) 

RETURN 
C 

998 

999 CALL ZZBCMB (5HYOUTI ♦ NERROR ) 

END 


NERROR=3 
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SUBROUTINE YPART (NUTA ,V,LV,KV,NUT1 ) 

DIMENSION V 1 1 > ,LV(1 ),MHEAD(10) ,MPHEAD(10) ,M2HEAD( 1C) 

DATA LIT, NOT/S, 6/ 

C 

C REPARTITION SPARSE MATRIX A BY ROWS. 

C USED TO REPARTITION SPARSE MATRIX A WHICH WAS FORMED IN A PROGRAM 
C HAVING A DIFFERENT DIMENSION KV. 

C THE MAXIMUM ALLOWABLE PARTITION SIZE IS KV/4 . 

C CALLS FORMA SUBROUTINES YIN ,YINI ,YOUT ,YOUTI t ZZBOMb. 

C DEVELOPED FY R A PHILIPPUS. NOVFMPFR 1968. 

C LAST REVISION BY WA BENFIELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA * LOGICAL NUMBER OF UTILITY TAPE C\' WHICH MATRIX A IS STORED. 

C V =■ VECTOR WOPK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NEPRCR EXPLANATION 

C 1 = SIZF (IMITATION FXCEEDED. 

c 2 = size Limitation exceeded. 
c 

LAE=0 

NREOG 

nnzpt=c 

DC 5 1-1,10 
5 M2HE AD ( I ) - 0 
MAXP=KV/4 

IF (MAXP.EC .0) M AXP=1 
LPPS=1 
REWIND NUTA 
REWIND NUT1 
C 

CALL YINI (NUTA ,MHEAD, 1,1C) 

KVCHK - MFTAD(6) 

IE (KVCHK. EQ.KV) RETURN 
NNZA = MHFAD(A) 

IF (NNZA. LE. 1 ) RETURN 
10 CALL YINI (NUT A *MPH E AD * 1 » 10 ) 

NNZP = MPHEADd ) 

NERR0R=1 

IF (NNZP.GT.KV) GO TO 999 
L AF = LAE +NNZP 
IF (NN2F.GT.C) GO TO 20 
CALL YINI (NUT A»LV»KV,KV) 

CALL YIN (NUT A,V»KV,KV) 

20 IF (LAE.GT.KV) GO TO 40 
30 LPRE=LPRS— l+NNZP 

NERR0R=2 

IF (LPRE.GT.KV) GO TO 949 
LAE=LP°e 

CALL YINI ( NUT A,Li/,LPRS,LPRE) 

CALL YIN (NUTA,V,LPRS,LPRE) 

NNZPT=NNZPT +NNZP 
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LPRS=LPPE+] 

IP (NMZPT .LT .NNZ A ) GO TO 10 

40 IF (LAE.GT.KV) L AE-LA E-NNZP 
L=MAXP 
LPWS = 1 
LPWE=L 

IF CLAF.LF.LPWF) GO TO 80 
50 IF (LPWE. FQ. LAE ) GO TO 80 

IF CLVt LPWE )/lCOOOO .LT. LVCLPWE+I ) /1C0000 ) GO TO 80 
IF (MAXP.LT .2) GO TO 80 

DC 60 I=2,MAXP 
L=L— 1 

IF ( LV C LPWE ) /I COO 00 .EC. LVIU/100000) GO TO 60 
L PWE=L 
GO TO eo 
60 CONTINUE 

80 IF ( LPWF.GT.LAE ) LPWE = LAE 
L-LPWF 

NNZPW-LPWE-LPWS+I 
M2HEADC1) = NNZPW 
M2HEAD ( 2 ) = LV(LPWS) 

M2HEADO) = LVCLPWE ) 

CALL YOUTI (NUTI *M2HEA0 ♦ 1 , 10 1 
NRFC=NREC+1 

CALL YOUTI (NUT! t LV» LPWS tLPWE ) 

CALL YOUT (NUTl,VtLPWS,LPWE» 

IF ( LPWE. EC. LAE . AND. NNZPT. EC? .NNZA ) GO TO 105 
IF (LAF.LT.LPWF+MAYP) GO TO 85 
LPWS=LPWE+1 
LPWE=LPWS— 1 +MAXP 
L=LPWF 
GO TO 50 
85 MOVE=LAE-LPWE 
LAE=MOVE 

IF (MOVE.EO.O) GO TO 1C 2 

DO 100 1=1 »MOVE 
L =L + 1 
V(I)=V(L) 

100 LV(I)=LV(L) 

102 LPP S = LAE + ! 

IF INNZPT.LT .NNZ A ) GO TO 3C 
L PWS= 1 
LPWE=LAE 
^0 TO 80 
105 REWIND NUT1 
REWIND NUTA 
MHFAD ( 3 ) = NREC 
MHEADC6) = KV 

CALL YOUTI (NUT A *MHE AD* 1 * 1C ) 
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DC 110 I=l t NREC 

CALL YINI (NUT1,MPHEAD»1»10) 
CALI YOUTI (NUT A tMPHE AD » 1 * 10 ) 
NNZP = MPHEADd ) 

CALL YINI (NUT1 ,LV, 1 ,NNZP) 

CALL YIN (NUT1,V,1»NNZP> 

CALL YOUTI (NUTA*LV*1.NNZP! 

CALL YOUT (NUTA,V,1,NNZP) 

110 CONTINUE 
C 

RETURN 

C 

999 CALL ZZBOMB (5HYP ART ,NERRCR) 
END 
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SUBROUTINE YPUNCH (NUTA « ANAME , V, LV, KV ) 

DIMENSION VII) ,LV(1),W(4) ,MHEADI10) 

DATA NIT ,NOT/5 ,6/ 

C 

C PUNCH SPARSE MATRIX A ON CARDS IN FORMA COMPATIBLE FORMAT FOR 
C SUBROUTINES RFAD*YRE AD* JREAD. 

C CAUTION - ELEMFNT LOCATIONS (LV) SHOULD BE IN INCREASING ORDER. 

C CALLS FORMA SUBROUTINES YIN *YINI *ZZBOMB. 

C DEVELOPFD PY R A PH1LIPPUS. JANUARY 1969. 

C LAST REVISION BY WA BFNF I ELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 
C ANAME = MATRIX IDENTIFICATION. (A6 FORMAT) 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V*LV IN CALLING PROGRAM. 

C 

C NERROR EXPLANATION 

Cl- SIZE LIMITATION EXCEEDED. 

C 

400! FORMAT ( A6,I4 V 15 , A6 ) 

4002 FORMAT <?15,4E17.8) 

4003 FORMAT ( IOHOOOOOOOOOO ) 

C 

REWIND NUTA 

CALL YINI (NUTA,MHEaD,1,I0) 

NR A = MHEAD(l) 

NCA = MHF AD ( 2 ) 

NPAPTA = MHEADm 
NNZA = MHEAD (4) 

ISHAPE = MHEAD ( 7 ) 

PUNCH 40C1* ANAME *NRA *NCA * ISHAPE 
TF (NNZA. EC;. 0) GO TO 40 
LPS= 1 
I FL AG=0 

C 

DO 38 M-l * NPAPTA 

CALL V1N1 (NUT A »MHE AD* 1*10) 

NNZP = MHFAD(l) 

LFELP = MHFAD(2) 

LLELP = MHE AD ( 3 ) 

IF (NNZP.GT.C) GO TO 2 
CALL YINI (NUTA *M, HEAD, 1 *2) 

CALL YIN (NUTA, V,l,2) 

GO TO ?.P 

2 LPE-LPS— 1 +NNZP 

NERR0R=1 

IF (LPE.GT.KV) GO TO 999 
CALL YINI ( NUT A*LV,LPS»LPF ) 

CALL YTN (NUTA,V,LPS,LPE ) 

r 

DO 3E I=LPS,LPF 
1 A=LV( I )/)CCCCC 
JA=LV ( I )-1 00000*1 A 
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IF (I.EC.l .AND. M.EQ.l) GC TO 2C 
K = JA-JS+1 

IF (IA.NF.1S .OR. K.LE.O .OR. K.GT.4) GO TO 5 
W(K)=VI I) 

IF (I.LT.LPE .OR. M.LT.NPARTA) GO TO 35 
IFLAG=I 
5 N J=4 

IF ( (JS+3) .GT.NCA ) NJ=NCA-JS+1 

IF (JA.GT.NCA) NJ-4 

PUNCH 4C0? t IS,J$*(W( J),J=1,NJ) 

IF (NNZP.EQ.O) GO TO 38 
IF (IFLAG.EO.I) GC TO 35 
20 I S=I A 
JS=JA 
W(1)=V(I) 

DO 30 K -? »4 
30 W(K)=0. 

IF (I.LT.LPE .OR. M.LT.NPARTA) GO TO 35 
I FLAG=1 
GO TO 5 
35 CONTINUE 
38 CONTINUE 

40 PUNCH 4003 
RETURN 

999 CALL ZZEOMB (6HYPUNCH ,NEERCR ) 

END 



onooooor»noooonnnnoor>or>oo<.fOoooooooonnnor>nooor>oooo 


YREAD 


1/ 5 


SUBROUTINE YREAD (NUTA , V » L V, KV ,NUT 1) 

DIMENSION vm,LV(1 )»X(4),IP.EMRK(9) ,MHE AD ( 1 0 ) 

DATA NIT, NOT/5 ,6/ 

READ SPARSE MATRIX EROM CARDS OR T APE AND PRINT IT. WRITE MATRIX ON 
TAPE IF SO INDICATED (EY NWTAPE .NF. 0 IN COLUMNS 79-BO). 

CALLS FORMA SUBROUTINES INTAPE, LT APE , PAGE HD , YIN ,YIN1 ,YCUT , 

YOUTI , YPART , YRTAPE , YWR ITE , YWTAPE ,2ZB0MB . 
DEVELOPED PY R A PHILIPPUS. NOVEMBER 1968. 

LAST REVISION BY WA BENFIELD. MARCH 1976. 


**** CARD INPUT **** 


FIRST CARD 


MIDDLE CARDS 


LAST CARO 


ANAME ,NROWS,NCOLS WITH A6,I4,I5 FORMAT. 

REMARKS IN COLUMNS 16-69. 

IF THIS IS THE LOWER (OR UPPER) HALF OF A SYMMETRIC 
MATRIX, THE WORD LOWER (OR UPPER) MUST APPEAR IN 
COLUMNS 16-2C. 

$ IN COLUMN 72 FOR NWTAPE INITIALIZATION. 

CONTRL IN 73-78. CONTRL MAY BE BLANK, OR THE WORDS 
REWIND OR LIST, OR (WHEN S IN 72) THE NWTAPE TAPE-ID 
(EG T1234). 

NWTAPE (LOGICAL TAPE NUMBER) IN COL 79-80 (EG 12). 
DATA WITH FORMAT (?15, 4EI7.0). 

1— ST 15 I S THE ROW NUMBER. 

2- ND 15 IS THE COL NUMBER OF THE NEXT E17.0 FIELD. 
NEXT 4E17.C ARE ELEMENTS OF THE MATRIX. 

TEN ZEROS IN COLUMNS 1-10. 


**** TAPE INPUT **** 

ONE CARD - ANAME, LOC(ZERO CR MINUS THE LOCATION NUMBER ON NRTAPF ) , 

NPTAPE (NUMBEP OF READ TAPE, IF MINUS NO PRINTOUT), ARUNNO 
PCNTRL(PLANK, REWIND, OR LIST) WITH FORMAT ( A6 , 14, 1 5 , A6, A6 ) 

- REMARKS IN 2 e- 69 . 

- $, CONTRL, NWTAPE SAME aS FIRST CARD L' k, DER CARD INPUT. 


SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA - LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 

V = VFCTOR WORK SPACE. 

IV = VFCTOR WORK SPACE. 

KV = DIMENSION S I ZF OF V,LV IN CALLING PROGRAM. 

NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

NEPROR EXPLANATION 

1 = ROW OR COLUMN VALUE OF CLEMENT FXCFEDS MATRIX SIZE. 

2 = DATA ON CARD PAST MATRIX COLUMN SIZE. 

3 = MATRIX IS DESIGNATED LOWER HALF SYMMETRIC BUT NON-ZEROES EXIST 

ONLY IN ThE UPPER HALF. 

4 = MATRIX IS DESIGNATED UPPER HALE SYMMETRIC BUT NON-ZEROES EXIST 

ONLY IN THE LOWFP HALF. 

5 = LOCATION ON TAPF PAST END OF TAPE MARK. 

6 = LOCATION ON TAPE PAST END OF TAPF MARK » 

1001 FORMAT (A6, 1A, 15 ,9A6, 2XA1,A6,I2) 

100? FORMAT C?H,4E17.0) 

2001 FORMAT ( //26H CARO INPUT SPARSE MATRIX A6 , ?X,IH(, 14, 2H X,I4,2h ), 
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* 

2002 FORMAT 

* 

2003 FORMAT 

2004 FORMAT 

2005 FORMAT 

* 

* 

2006 FORMAT 


2X 9A6,2X A1,A6,I4//) 

(//?6H CARO INPUT SPARSF MATRIX A6 * 2X V 1H(,I4,2H X,I4,2H ) 
3X, 9HC0NTINUFD //) 

(// 1 XA6 » 1 4, 1 5 » 5X 9A6,2X Al,A6,I4) 

(1 X 215 »4f 17.8 ) 

(14H0FND OF YREAD. 

/31H NUMB FR C'F NON-ZERO ELEMENTS * 15 
/24H NUMBER OF PARTITIONS = 13) 

(25H0SIZE OF MATRIX READ IS (14, 2H XI4,2H )) 


READ IN HEADER CARD. 

NUT=NUT1 

READ (NIT, 1001) ANAME ,N1 ,N2, IREMRK , IZ 1 , IZ2 *NWTAPE 
CALL PAGEHD 


IF (N1.LE.0) GO TO 200 


CARD READING SECTION. 

L0W=0 

LUP-0 

REWIND NUT 
NNZ A=0 
NR EC = 0 
LAE=0 
NR A=N1 
NCA=N2 

WRITE (NOT, 2001) ‘ ANAM E,NRA ,NC A , IREMRK ,IZ1 , IZ2 ,NWTAPE 
NlINE=0 

11C READ (NIT, 1002) I ,JS, X 

IF (I.EO.O .AND. JS.EO.O) GO TO 132 

NERR0R=1 

IF (I.LF.O . OR • l.GT.NRA .OR. JS.LE.O .OR. JS.CT.NCA) GO TO 990 
JE = JS<3 

IF (JF.LE.NCA) GO TO 115 
JX=NC A— JS+2 

NERRCR=2 

DO 112 J=JX,4 

IF <X(J) .NE. 0.) GO TO 990 
112 CONTINUE 
JF = NC A 
115 N = 0 

DO 120 J = JS , JE 
M - N+l 

IF (X(N).FO.O.) GO TO 120 
IF (I.CT.J) LOW= 1 
IF (I. | T. J) LUP= 1 

IF (IREMRK (1 ). F 0.5HLUWER .AND. I.LT.J) GO TO 120 
IF (3REMPM 1 ).EQ. 5HUPPER .AND. I.CT.J? GO TO 120 
LAE = L AE +1 

IF (LAE .LF.KV/4) GO TO 118 
LAE = LAE-1 
MHEAD(l) = LAE 

CALL vpUTT (NUT 1 ,MHE AD, 1 , I ) 

NRFC“NR P C *1 

CALL YOUT! (NUT1 ,LV , 1 ,L AE ) 
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CALL YOUT (NUT1 »V» 1 »LAF) 

L AE = 1 

118 NNZA=NNZA+1 
V ( LAF ) =X ( N) 

LV(LAF )=100000*I+J 
l?o continue; 

NLINF - NL1NE+1 
IF (NLINE.LE.47) GO TO 1*5 
CALL PAGFHD 

’■'RITE (NOT *2002 ) AN AM E »NR A »NCA 
Nl INF = 1 

125 WRITE i NOT ,2004 ) I,JS,X 
GO TO 110 

132 IF (LAE. EC. 0) GO TO 135 
MHEADH ) = LAF 

CALL YOUT I (NUT 1 ,MHE AD ,1,1) 

NRFC=NR FC+1 

CALL YCUTT (NUT1 ,LV, 1 *LAE > 

CALL YOUT (NUT1 ,V,1 ,LAE) 

135 REWIND NUT 
REWIND NUT A 
IASHAT= JRFMRKP ■ 

IF ( I ASHAP.NE.5HWHCLE .AND. IASHAP.NE .5 HUP PER .AND. 

♦ IAFHAP.NF.5HL0WER .AND. IAF^t AP.NE.4HDI AG) 1ASHAP=5HWH0LE 
IF (NNZA.EO.O) GO TO 137 

NERR0R=3 

IF (IREMRM1 ) • EG * 5HL0WER .AND. LOW.FO.C .AND. LUP.EC.Il GO TP 999 

NFPR0R=4 

IF (IRFMRKd ).F0. 5HUP PER .AND. LUP.EO.O .AND. LOW. EC. 1) GO 7t‘ 999 

137 MHE AO ( 1 ) = NRA 
MHE AD (2) « NCA 
MHE AD ( 3 1 = NREC 
MHFAD(4> = NNZA 
MHE AD ( 5 I = 0 
MHF AD ( 6 ) = C 

MHE AD ( 7 ) = IASHAP 
MHE AD ( 8 ) * 0 
MHF AD ( 9 ) = 0 

CALL YOUT I I NUT A *MHE AD, 1 , 10 ) 

DC 139 1=1,10 
139 MHF AD ( I ) - 0 

IF (NNZA.GT.O) GO TO 138 
CALL YOUT I (NUTA,MHEAD,1,1C) 

CALL YOUT I (NUT A ,MHE AD , 1. » ? 1 
CALI YOUT 1 (NUT A ,MHE AD, 1 »2 ) 

GO TO 3C0 

138 DC 1*0 i=l, NREC 

CALL YIN I (Niri ,MHEAD,1,1) 

NNZP = MHE AD (1) 

CALL YINI (NUT 1 ,LV , 1 ,NNZP ) 

CALL YIN (NUT l ,V,1 ,NNZP ) 

MHE AD ( ? ) = L V( 1 ) 

MHF AD ( 3 ) = LV(NNZF) 

CALL YOUT I (NUTA ,MHE AD, 1 , 10 ' 

C ALL YOUT 1 (NUT A ,LV , 1 »NNZP \ 
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140 CALL YOUT INI, fA,V,l ,NNZF « 

CALL YPART (M'TA,V,L V,KV,NUT1 ) 

GC TO 300 

TAPE POSITIONING AND READING SECTION, 

200 WRITE <N0T,?CC3) *NAME ,N1 ,N2 * I REMRK , IZ1 »IZ 2, NWTAPE 
NRTAPE = IAFS(N2) 

Tr HREMRM2) .EQ. 6HREW1ND) REWIND NRTAPE 
IF IIREMRM2) .EQ.AHLIST) CALL LTAPE I NR T APE ) 

IF CN1.FQ.0) GC TO 250 
POSITION NRTAOf. 

READ (NRTAPE) T ID,LN , IEOTCK 
NUM = LN+N1 

NERR0R=5 

IF (NUM) ?C5,??0,2?5 
205 IF I IECTCK . EC.3HECTI GO TO 900 
READ ir’RTAP E) 

NUM = -NUN-1 

IF (NUM.EQ.O) GO TO 25C 

NERR0R=6 

DO 210 L=1,NUM 

RF AD (NCTAPE) TIC.LN* .ofCK 
IF (IE0TCK.EC.3HE0T) GC TO 900 
210 READ (N*?TAPE ) 

GC TO 250 

220 BACKSPACE NRTAPE 

GO TO 250 

NOTE-. .THE FCLICWJNG SECTION WAS DESIGNED PRIMARILY TC BE USED ON A 
DISK. IT WILL WORK ON A TAPE BUT IT WILL NOT BE AS EFFICIENT 
********* * 

225 REWIND NRTAPE 
NUM = ( — N 1 - 1 )*2 
IF (NUM.EC.ei GO TO 250 
DO 230 L - 1 » NUM 
230 READ (NRTAPE) 

********** 

250 CALL YP.TAPE ( 1REMRK ( 1 ) ,ANAME , NUTA, V ,LV,KV, NRTAPE *NUT1 ) 

REWIND NUT* 

CALL V INI (NUT A ,MHE AD, 1 , 10 ) 

NR A = MHEAD(J) 

MCA - MCE AD (2) 

NPART = MHEAD( 3 ) 

NNZ A = MHEAD(A) 

WRITF ( NOT ,2CCfc ) KRA.NCA 

IF (N2.CT.0) CALL YWRITE ( NUTA ,AN AME , V, LV , KV) 

NWTAPE INITIALIZING, WRITING, AND LISTING SECTION. 

300 IE (NWTAPF.LF.O) GC TO 350 

IF (IZ1.F0.IHt ) r f- 1 L INTAPE (NWTAPE, IT. 2 ) 

IF (IZ?.FC.6HRFWIN,'KEWIND NWTAPF 
CALL VWTAPF (NUTA ,ANAME ,V, LV ,K V, NWTAPE ) 

IF UZ2.EC..AHLIST) CALL LTA? (NWTAPF) 


350 REWIND NUTA 
CALL YIK'I 


(NUT A ,MHE AD* 1 , 1C 
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NR A = MHEAD(l) 

NCA = MHEAPI2) 

NPART = MHEADfS) 

NNZ A = MHBAD(A) 

WRITE (N0T*2C05 1 NNZA.NPART 
RETURN 

900 CALL LTAPE <NRTAPE> 

CALL ZZPOMB I5HYREAD t NERROR) 
990 WRITE (NOT,! PM I,JS*X 
9°9 CALL ZZBOMB (5HYREAD *NERR0R ) 
END 
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SUBROUTINE YKLVAD ( ALPHA, NUTA, I VEC • JVEC»NUTZ » V,LV,KV,NUTl ,NUT2 * 

* NUT3,NUT4) 

DIMENSION V(1),LV(1), IVEC(1),JVEC(1 ),MHEAD(10) 

OATA NIT *NOT/S *6/ 

REARRANGE ROWS AND COLUMNS OF ALPHA * SPARSE MATRIX A AND ADD TO 
SPARSE MATRIX Z. (ALPHA * A ♦ Z = ZI. 

CALLS FORMA SUBROUTINES XLORD * YA ABB ,YIN ,YINI ,YLC>RD ,VNCZER, 

'OUT * YOU TI »YPART » YSYMLH, YSYMUH*ZZBOMB. 
DEVELOPED PV R A PHILIR^US. DECEMBER 1969. 

LAST REVISION PY RL WOHlRN FOR NASA. MAY 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

ALPHA = SCALAR THAT MULTIPLIES MATRIX A. 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 
IVEC = VECTOR. (SIZE = NR A) 

TVECU) = ROW POSITION OF A(PCW I) IN Z. 

IF IVFC(I) IS PLUS , Z = Z(POW IVEC ( I ) ) * ALPHA * A(ROW I) 

IF IVFC(I) IS MINUS, Z = Z ( ROW IVEC(I)) - ALPHA * A(ROW I) 

IE ( IVEC ( I ) IS ZERO , A (ROW I) IS OMITTED IN Z. 

JVEC = VFCTCR. (SIZE = NCA) 

JVEC(J) = COLUMN POSITION OF A (COL J) IN Z. 

IF JVEC(J) IS PLUS , Z = Z ( COL JVEC(J)) ♦ ALPHA * A(COL J) 

IF JVFC(J) IS MINUS, Z = Z(COL JVEC(J)) - ALPHA * A(COL J) 

IF JVEC(J) IS ZERO , A (COL J) IS OMITTED IN Z. 

NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZF OF V,LV IN CALLING PROGRAM. 

NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

NUT 2 = LOGICAL NUMBER OF UTILITY TAPF. 

NUT3 = LOGICAL NUMBFR OF UTILITY TAPE. 

NUT4 = LOGICAL NUMBER OF UTILITY TAPE. 

NFRPO*> EXPLANATION 

1 = ROW LOCATION OUT SIDE MATRIX Z. 

2 = COLUMN LOCATION OUTSIDE MATRIX Z. 

CALL YPART (NUT A ,V,L V,KV,NUTI ) 

C GET (A) MFADFR INFORMATION. 

REWIND NUTA 
REWIND NUT1 

CALL VTNI (NUTa,MHEAD,1,10) 

NR A = MHFAOm 
NCA = MHE AD ( 2 ) 

NPARTA - MHE AD (3) 

NNZ A = MHEAD (4 ) 

IE (NNZA.EO.O) PF* " N 
ISHAP - MHF AD ( 7 ) 

C GFT (Z) HEADER INFORMATION. 

REWIND NUTZ 

CALL YINI (NUTZ,MHEAL ’,?) 

NRZ = MHE AD (1) 

NCZ = MHF AD ( ? ) 
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jC CHECK SIZES. 

NERR0R=1 

DC 2 1=1, NPA 

IF 1 1 APS ( IVEC( I!) .GT.NR2 I GO TO 999 

2 CONTINUE 

NERR0R=2 

DO 3 1=1, NCA 

IF <IABS( JVECU ) i .GT.NCZ) GO TO 990 

3 CONTINUE 
C 

IF (ISHAP.NE.4HDIAG) GO TO 10 
ISHAP=5HWHCLF 
IF (NRA.NE.NCA) CO TO 10 
C 

DO 5 1=1, NRA 

IF CIABSaVECm ).NE.IABS< JVECCim GO TO 10 
5 CONTINUE 

c 

I$HAP=4HDIAG 
10 MHEA0I5) = 0 
MHEAO (6 ) = 0 

CALL YOUTI (NUT1 ,MHE AD, 1 , 10 ) 

C 

C BLOW-UP (A) TO (Z) SIZE. 

DC ICC i=i,nparta 

CALL YINI (NUTA,MHEAD,',10) 

NNZPA = MHEAD ( 1 ) 

CALL YINI (NUTA ,LV,1, NNZPA ) 

CALL YIN (NUTA,V,1 , NNZPA) 

C 

DO 50 J = 1 , NNZPA 
IA=LV(J)/1PCC00 
JA=LY ( J )— 100000*1 A 
IF (IVFC(IA)) 15,25,35 
15 IF (JVEC(JA)) 20,25,30 
20 LV(J)=-1CC0CC*IVEC(IA)-JVEC( JA) 

GO TO 5C 
25 V(J)=0. 

GO TO 5C 
30 V(J)=-V(J) 

LV( J) =-100000*1 VEC( I A)+JVEC( JA) 

SO TO 5C 

35 ir (JVECIJA): 40,25,45 
40 V(J)=-V(J) 

L V ( J ) = 1 00000*IVFC ( I A )— JVEC ( JA ) 

GC TO 5C 

45 LV( J)=ICOOOO*IVEC (IA) ♦JVEC(JA) 

50 CONTINUE 
C 

IF ( I SNAP • EU. 5HWH0L E .OF. IS HAP. E 0. 4HDI AG ) GO TO 9C 
IF (ISHAP.E0.5HL0WER ) GC TC 70 
c 

DC 60 K=1 , NNZPA 
1A=LV!K)/1CCC0C 
JA = LV ( K )— 1 OCOCO* I A 
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IF (IA.GT.JA) LV(K)=1C00G0*JA+IA 
60 CONTINUE 
*C 

GO TO 9C 
C 

70 DC 80 K=I,NNZPA 
lA=LV(K)/lCCrOO 
JA=LV(K ) — 100000*1 A 
IF (IA.LT.JA) LV(K)=100000*JA+IA 
80 CONTINUE 
C 

90 MHEADC2 ) = LVC1) 

MHE AD ( 3 ) = LV(NNZPA) 

CALL YOUTI (MJT1 -MHEAD,1,1C) 
CALL YOUTI (NUT1 .LV.l.NNZFA) 

100 CALL YCUT (NUT1 ,V, 1 ,NNZPA) 

C 

CALL YNOZEP (NUT1 ,V,LV,KV,NUT2 ) 
REWIND NUTZ 
REWIND NUTZ 

CALL YINI CNUTZ.LV. 1 ,10) 

CALL YOUTI (NUT2.LV, 1,10) 
NPARTZ-LV (3 ) 


TRANSFER HFIGINAL (Z> FROM NUTZ TO NUT2. 
DO 110 J= 1 ,NPAR T Z 
CALL YINI (NUTZ ,LV, 1,10) 

(NUT2.LV, I, IC) 


CALL YOUTI 
NNZ=L V ( 1 ) 
105 CALL YINI 
CALL YIN 


(NUTZ ,LV , 1 ,NNZ ) 
(NUTZ ,V, 1 »NN2 ) 


CALL YOUTI (NUT2.LV. l.NNZ) 
CALL YOUT (NUT 2 ,V , I ,NNZ ) 
110 CONTINUE 


ADD ALPHA*? I OWN-UP (A) TO ORIGINAL (Z). 

CALL YAAEB ( ALPHA, NUT1 , 1 . ,NUT2 ,NUT2 ,V, LV , KV, NUT3 ,NUT4 ) 
RETUFN 

°99 CALL ZZ50MB (6HYP FVAD ,ME RROR ) 

END 
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C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


SUBROUTINE YRTAPE ( I A RUNG* IANAMF »NUTA »V » LV t KV » NTA PE »NUT1 ) 
DIMENSION VM)»LV(!)t MCHECKC 2 ) *MHFAD( 10 ) 

DATA NIT*N0T/5,6/ 

RFAD SPARSF MATRIX A FROM TAPE ( NTAPE ) BY IDENTIFICATION OF IARUNC 
AND IANAMF AND STPRF IT ON UTILITY TAPE (NUTA ) • 

CALLS FORMA SUBROUTINES LTAPE ,YIN ,YINI ,YOUT ,YCUTI .ZZBOMB. 
DEVELOPED FY R A PHILIPPUS. NOVEMBER 1968. 

LAST REVISION BY KA BENFIELD FOR NASA. MAY 1976. 

SUBROUTINE ARGUMENTS CALL INPUT) 

IARUNO - RUN NUMBER OF MATRIX A. (A6 FORMAT) 

IANAMF - MATRIX IDENTIFICATION. CA6 FORMAT) 

NUTA = LOGICAL NUMBEP OF UTILITY TAPE ON WHICH MATRIX A IS STORED. 
V = VFCTOP WORK SPACF. 

LV = VFCTC'R WORK SPACE. 

KV - DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

NTAPE = LOGICAL NUMBER OF TAPE FROM WHICH MATRIX A IS TO BE READ. 

NUT1 = LOGICAL NUMBEF OF UTILITY TAPE. 

NEPROR EXPLANATION 

1 = INCORRECT MATRIX TYPE. 

2 = DIMENSION SIZE EXCEEDED (KV). 

3 * IANAME AND IARUNO NOT FOUND ON NRTAPE. 

4 = DIMENSION SIZE EXCEEDED (KV). 

3001 FORMAT (3 OH 1 YRTAPE CANNOT FIND RUNNO = A6/ 

* 22X8HANA.KE = A6/17X I3H PARTITION NO. I5/30X6H~ ) 


NPF = 1 
NREC=C 
NNZ A=0 
NTIMF=0 
REWIND NUTA 
REWIND NUT1 

SEARCH TAPF FOR CORRFCT HEADING. 

5 READ (NTAPE) TAPE ID ,LN , I EOTCK t ITRUNO, ITNAME , NR A»NCA, DATE * ITYP ,NNZP 
*, NP , NPT t ( MCHECK (I ) * 1=7 »2),ISHAP 

IF (I5HAP.FC.0 .AND. NRA.NE.NCA) I SHAP=5HWHCLE 
IF ( ITRUNO. FQ. IARUNO .AND. ITNAME. EQ. IANAME ) GO TO 1C 
IF (IFPTCK.FQ.3HE0T) GO TO 20 
READ (NTAPE) 

GO TO 5 

MATRIX HAS EE.FN FOUND. 

10 NERRCR= 1 

IF ( ITYP. NE ,5H$P ART .AND. ITYP ,NF .6HS PARSE ) GO TO 99C 
IF (ITYP.NE.BHSPART) GO TO 32 
IF (NP.FO.NPF) GO TO 1? 

RFAD (NIAPE) 

GO TO 5 

12 NERRCR=2 

IF (NNZP.GT.KV) GO TO 99C 
IF (MNZP.GT .0 .OR. NPT.EQ.l) GO TO 15 
READ (NTAPE) 

NPE =NPF + 1 
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IF (NPF.GT.NPT) GO TO 25 
GO TO 5 

15 RFAD (NTAPE) ( LV ( I) * V ( I ) ,1 =1 .NNZP ) 
MHFAD(l) = NNZP 

MHfc ADI 2 ) - LV(J) 

MHEADO) = LV(NNZP) 

DO 16 I =4* 1C 

16 MHFADII) = 0 

CALI. YCUT1 (NUT1 .MHEAD.I.1C) 

CALL YOUTI (NUTl.LV.l.NNZP) 

CALL YOUT (NUTl.V.l .NNZP) 

NTIME=0 

NRFC-NREC +1 

NPF=NPF+1 

NNZA=NNZA+NNZP 

IF (NPF.GT.NPT) GO TO 25 

GO TO 5 

C SEE IF MATRIX WAS FOUND. 

20 NTIME=NTIME+1 

IF (NTTMF.F0.2) GO TO 900 
REWIND TAPE 
GO TO 1 

25 MHEADO ) = NRA 
MFEAD ( 7 ) = NCA 
MHE AD ( ? ) = NR EC 
MHEADC4) = NNZ A 
MHE AD » 5 ) =• MCHECK (1 ) 

MHEADC6) = MCHECKI2) 

MHEAD(T) = ISHAP 
CALL YOUTI (NUTA.MHEAD.1.10) 
REWIND NUT1 
C 

DO 30 1=1 .NREC 

CALL YINI (NUT1 .MHEAD.l.lD 
CALL YINI (NUT1 .LV. 1 , MHE AD ( 1 ) ) 
CALL YIN (NUTl.V.l .MHEADCl)) 
CALL YOUTI (NUTA.MHEAD.1,10) 

CALL YOUTI (NOT A .LV, 1 ,MHEAD( 1 ) I 
30 CALL YOUT (NUT A .V.l ,MHEAD(1 ) ) 

C 

RETURN 

32 

IF (NNZP.GT.KV) GO TO 990 

RFAD CNTAPF) ( LV ( I) .V ( I ) .1 = 1 .NNZP) 

MHEADd ) = N' A 

MHF AD ( 2 ) = NCA 

MHEADO ) = NPF 

MHFAD(4) = NNZP 

MHE AD ( 5 ) = MCHECK !l) 

MHE AD ( 6 ) = MCHECK (2) 

MHE AD (7) = ISHAP 

CALL VOUTI (NUTA.MHt AD, 1,10) 

MHEADd ) = NNZP 

MHE AD (2 ) = LVd) 


NEkRCR=3 


NERR0R=4 
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MHEAD(3) = LV(NNZP) 

MHFaD(4) = 0 
MHEADI5) = 0 
MHEAD (6 ) = 0 
HHEADC7 ) = 0 

CALI YOIJTI (NUT A ,MHEAD» 1 1 10 ) 

CALL YOUTI (NUTA »LV» 1 »NNZP ) 

CALL YOUT (NUT A ,V « 1 ,NNZP ) 

RETURN 

900 WRITE (NOT, 3001) IARUNP, IANAME ,NPF 
990 CALL LTAPE (NTAPE) 

CALL ZZBOKB (6HYRTAPE *NERROR ) 

END 
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SUBROUTINE YRV1 <NUTI,N,NU,7,LV,KV,NUT1 ,NUT2,NUT~ ) 

DIMENSION V(l), LV(1) 

DIMENSION IH(IC) 

DATA IRtIL / 100, lOOOOOOOOCOCO / 

DATA IMULT / 899999999 / 

DATA SCALE / 9999999999. / 

C 

C GENERATE MATRIX OF RAYLEIGH VECTORS. 

C RANDOM NUMBERS BF TWEEN -1. AND +1.. 

C INPUT NON-ZERO COLUMNS ARE NOT CHANGED. 

C CALLS FORMA SUBROUTINES ... 

C YIN , YINI , YL OPD ,YOUT ,YOUTI , YPART , YTRANS , ZZBC'MB . 

C DEVELOPED BY RA PHILIPPUS. MARCH 1972. 

C LAST REVISION BY RL WOHLEN. JUNE 1975. 

C 

C SUBROUTINE ARGUMENTS fALL INPUT) 

C NUTZ = LOGICAL NUMPEP OF UTILITY TAPE OF MATRIX OF INITIAL RAYLEIGH 

C VFCTORS. NON-ZERO COLUMNS INPUT TO THIS SUBROUTINE ARE NOT 

C CHANGED. 

C N = NUMBER OF ROWS OF RAYLEIGH VECTORS MATRIX. 

C NU = NUMBFR OF CCLS OF RAYLEIGH VECTORS MATRIX. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION OF V,*V IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPF. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

‘ NUT3 = LOGICAL NUMBER OF UTILITY TAPE. 

K VI A = KVM 

C TRANSFER ROWS ON NUTZ TO COLUMNS ON NUT3 . 

CALL YTRANS ( NUTZ .NUT 3 , V»LV, KV ,NUT1 ,NUT2) 

REWIND NUT3 
REWIND NUTZ 
WW11 49WV 

NPARTZ = (NNZZ-1 ) /KV1A+1 
C SET (Z) HEADER. 

IH ( 1 ) = N 
IH ( 2 ) = NU 
IF ( 3 ) = NPARTZ 
IH(4) = NNZ Z 
IH(5 ) = C 
IH(fc) = C 
IH ( 7 ) = 5 HWHOLE 
I H( 6 ) = 0 
IH19) = 0 
IH(10) = 0 

CALL YOUTI (NUTZ, IH, I ,10) 

LZ = 0 

IRN = IMULT**? 

C READ DATA THAT ENTER FD ON NUTZ BUT TRANSPOSED TO NUT 3. 

CALL YINI (NUT3,IH, 1 ,10) 

NPARTI = IK( 3 ) 

NR FAD = 0 
U— 0 49 T 3 
LIS = L IE +1 
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LI = LIS 
NNZPI = TH(A) 

IZ = 0 
IZP = C 
IYCUT = 0 
DO 10 I=1,KV1A 

10 veil = 0 . 

C J=CCLUMN NUMBFR. 

DO 99 J=1,NU 

IF (NPARTI.EO.C) GO TO AO 
20 IF (NRFAD.EQ.NPARTI .AND. LI. GE. LIE I GO TO AO 
IF (LI .LE .LIE) GO TO 30 
CALL YINJ (NUT3,IH,1,10) 

NNZPI = 1H(1) 

LIE = LIS-l+NNZPI 

CALL YINI (NUT3 ,LV, L IS »LIE ) 

CALL YIN (Nl«T3, V, LIS, LIE) 

NR FAD = NREAD+1 
LI = LIS 

30 IZ = L''(LI)/100000 
IF (IZ.GT.J) GO TO AC 
IF (IZ.LT.J) L I=L I+l 
IF (LT.GT.LIF) GO TO 20 
IF (IZ.LT.J) GO TO 30 
JZ = LV(LI)— 10CCCO*IZ 
IF (IZ.EO.IZP) GO TO 3A 
LZ PI = L2+1 
LZPN = LZ+N 
LZPE = LZ FN 

IF ( LZPN.GT »KV! A ) LZPN=KV1A 

L = 0 

DO 33 K=LZP1,LZPN 
L = L+l 

33 LV ( K ) = 1 00000*L + J 
3A IF (LZ+JZ.LE.KV1A) GO TO 38 
IH ( 1 ) = KV1A 

IH(?) = LV(I ) 

IH ( 3 ) = LV(KVIA) 

CALL YOUTT (NUTZ ,IH, 1 ,10) 

CALL YOUTI (NUTZ»LV,1 ,KV1A) 

CALL YCUT (NUTZ, V,1,KV1A) 

IYOUT = IYOUT+1 
DO 3? K-l.KVIA 
L = L 4 -) 

LV(K) = 100000*L+J 

37 V ( K ) = 0. 

LZ = 1-JZ 

38 V ( LZ+JZ ) = V(LI ) 

IZP = IZ 

LI = L I + 1 

IF (LI.GT.LIF) GO TO Ai> 

GO TO 30 

AO IF (IZP.NF.J) GO TO Bo 
A5 LZ = LZ+N 
GO TO 99 
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C PLACE RANDOM NUMBERS IN ONE COLUMN. 

80 DO 85 1-J f N 

IRN = IRN*IMULT 

JRN = ( IRN— ( IRN/I L)*I L I/1R 

IF t(JRN/2t*2 .EC. JRN ) JP.N=-JRN 

RNUM = JRN 

LZ = LZ+l 

VCLZ) = RNUM/SCALE 
LV(LZ) = 100000*1 +J 
IF (LZ.LT.KV1M GO TO 85 
IHC1) = LZ 
IHf ? ) = LVII) 

IHC3) = LV(L2) 

CALL YOUTI (NUTZ ,IH, 1,10) 

CALL YOUTI <NUTZ,LV,I,LZ) 

CALL YOUT (NUTZ »V ,1,LZ) 

IYCUT = IYOUT+I 
LZ = 0 
85 CONTINUE 
99 CONTINUE 

IF ( I YOUT .EQ.NPARTZ ) GO TO 1C9 
IHtl) = LZ 
IH ( ?) e LV(1) 

IH ( 3 ) = LV(LZ) 

CALL YOUTI (NUTZ »IH » 1,10) 

CALL YOUTI (NUTZ,LV,1,L2) 

CALL YCUT (NUTZ » V,l,LZ) 

109 CALL YLORD (NUTZ,V,LV,KV,NUT1 ,NUT2> 
RETURN 
END 
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SUBROUTINE YRVAD1 ( AL PHA * A ♦ I JVEC ,NUTZ ,NRA * V,LV,K V ,KA,NUT1 ,NUT2 « 

* NUT3,NUT4) 

DIMENSION V(1),LV(1),IJVEC(1 ), MHEAD(IO), A(KA*1 ) 

DATA NIT »NC'i /5 ,6/ 

DATA EPS / l.E-25 / 

REARRANGE ROWS AND COLUMNS OF ALPHA * MATRIX A AND ADD TO 
SPARSE MATRIX Z. (ALPHA * A ♦ Z = Z). 

MATRICES A, 2 ARE ASSUMFD SYMMETRIC. 

CALLS FORMA SUBROUTINES XLOPD , YA»VJB ,YIN »YINI » YLORD ,YNOZER, 

YCUT » YOUT I , YPART , YSYMLH, YSYMUH,ZieOMB. 
DEVELOPED BY R A PHILIPPUS. NOVEMBER 1972. 

LAST REVISION FY WA EENFIELD. MATCH '976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

ALPHA = SCALAR THAT MULTIPLIES MATRIX A. 

A = MATRIX A. 

I JVEC = VECTOF. (SIZF = NR A ) 

IOVEC ( I ) - ROW POSITION OF A(ROW I) IN Z. 

IF I JVEC ( I ) IS PLUS , Z = Z ( ROW IJVEC(I)) ♦ ALPHA * A(ROW I) 

IF I JVEC ( I ) IS MINUS, Z = ZJROW IJVEC(I)) - ALPHA * A(ROW I) 

IF ( I JVEC ( I ) IS Z E p O , A ( ROW I) IS OMITTED IN Z. 

NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 
NR A = NUMBER OF ROWS AND COLUMNS OF A. 

V = VFCTCR WOPK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

KA = ROW DIMENSION OF MATRIX A IN CALLING PROGRAM. 

NUT! = LOGICAL NUMBER OF UTILITY TAPE. 

NUT? = LOGICAL NUMBER OF UTILITY TAPE. 

NUT3 - LOGICAL NUMBER OF UTILITY TAPF. 

NUT4 = LOGICAL NUMBER OF UTILITY TAPE. 

NFRROR EXPLANATION 

1 = 2 MATRIX IS NOT SQUARE. 

2 = ROW OP COLUMN LOCATION OUTSIDF MATRIX Z. 

3 * DIMENSION SIZE EXCEEDED ( KV ) . 

REWIND NUT1 
REWIND NUTZ 

CALL YINI (NUTZ,MHEAD, 1,10) 

NR2 = MHEAD(l) 

NCZ = MHF AD ( 2 ) 

NERROR-1 

IF (NRZ.NE.NCZ) GO TO 999 
ISHAP - MHFAD(7) 

IF ( ISHAP .EQ.9HL OWEF .OR. ISHAP. EQ. SHUPPER ) GO TO 1 
CALL YZFRUH (NUT2 ,V,LV,KV,NUTl ,NU72 ) 

REWIND NUT1 
REWIND NUTZ 

CALL YJNI JNUT2 ,MHE AD, I , 1 C ) 

ISHAF = SHLCWfR 
MHt ADf 7) = ISHAP 


NE RROR"2 
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DO 2 1 = 1 * NR A 

IF CIARSCIJVECU n.GT.NRZ) GO TO 999 
2 CONTINUE 

J = 0 

NERR0R=3 

DC 50 I A= 1 ,NR A 
DO 50 JA= 1* NR A 

IF (1A.GT.JA .AND. ISHAP.EC. 5HUPPER ) GO TO 50 
IF UA.LT.JA .AND. I SHAP .EC 1 • 5HL0WE R ) GO TO 50 
IF <ABS!A(IA,JA) I.LT.EPS) GO TO 50 
J = J+l 

IF ( J.GT.KV) GO TO 999 
IF (IJVECUA )) 1 5 *25 , 35 
15 IF (IJVEC(JA)) 20,25,30 
20 LVC JI = —I 00000*1 JVEC ( IA I—I JVEC ( J A) 

V(J) = A(IA,J«> 

GO TO 50 
25 J = J-l 
GO TO 50 

30 VIJ» = —A ( I A , JA ) 

L VC J ) = -1C0CCC* I JVEC ( IAI+IJVECC JA) 

GO TO 50 

35 IF CIJVECIJA)) 40,25,45 
40 V(J) = —A ( I A , J A ) 

LVCJ) = 100000*1 JVEC C lA )-?. JVEC CJA ) 

GO TO 50 

45 LV( J ) = 1 CC 000*1 JVEC l IA) U JVEC (JA ) 

V(J) = A ( IA , JA ) 

50 CONTINUE 

IF (ISHAP.F0.5HL0WER) GO TO 70 

DC ■:<' K = 1,J 

IA - LV(K )/l 00000 
JA = LV(K )-10C0C0*IA 
IF UA.GT.JA) LV!K)=100000*JA+IA 
60 CONTINUE 

GO TO 90 

70 DO EO K = 1 , J 

IA = LV (K I/10000C 
JA = LVCK >-100000*1 A 
IF (IA.LT.JA) LV(K)=100C0C*JA+IA 
80 CONTINUE 

90 MHE ADC?) = 1 
MHE AD (4 ) r J 
MHEADC5) = 5H0RDER 
MHE ADC 6) r v V 

IF C J.C-T.KV/ I MHFADC 6 ) = 0 
CALL XL OF C V,LV,1,J) 

CALL YC'ni (NUT1 ,MHE AD,1 ,10 ) 

MHE AD 5 * ! = J 



YRVAD1 


3/ 3 


MHFAD( 2 ) = L V ( 1 > 

MHE AC ( 3 ) = LV( J ) 

MHEADI4) = C 
MHEAD(5> = 0 
MHFAD (6 ) = 0 
MHFAD(7) = 0 

CALL YOUTI (NUTl , MHFAD, I ,10) 

CALL YOUTI (NUTl ,LV, 1, J) 

CALL YOUT (NUTl ,V, 1 , J ) 

C 

REWIND NUTZ 
REWIND NUT2 

CALL YIMI (NUTZ ,LV, 1 ,10) 

CALL YPUTI (NUT? *LV ,1,10) 
wPARTZ- LV (3 ) 

C 

DO 110 J=1,NPARTZ 

CALL YINI CNUT2 ,LV,1,10) 

CALL YOUTI (NUT? ,LV, 1,10) 

NNI=LV(1 ) 

IF (NNZ.GT-C) GO TO 1 05 
CALL YINI (NUT2 ,LV, 1,2) 

CALL YOUTI (NUT2 ,LV, 1,2) 

CALL YIN (NUTZ ,V,1 ,2) 

CALL YOUT (NUT2 ,V, 1 , 2 ) 

GO TO nr 

1C5 CALL YINI (NUTZ ,LV, 1 ,NN2 ) 

C»,lL YIN (NUTZ,V,1 ,NN7.) 

CALL YOUTI (NUT? ,1 V, 1 ,NNZ ) 

CALL YOUT (NUT 2 ,V,1 ,NNZ) 

110 CONTINUE 
C 

CALL YAAEE ( AL PH A, NUTl , 1 . ,NUT?,NU‘ \£ , V, LVt KV, NUT3 ,NUTa ; 
RETURN 
C 

999 CALL ZZBOMB (6HYRVAD1 ,NERROR ) 

END 



ooooooooooo 


YRVAD2-- 1/ 7 


SUBRCUT INF YRVAD 2 (NUTA, NUTZ ,NRZ ,W ,KW,V ,LV ,KV *NU7 1 ,NUT2 ,NUT3 ) 
DIMENSION V(1),LV(1),IU(16),IL(16),MHEAD(1C),MPHEAD<10>,M2HEAD<10) 
D IMENS I C’N W(KWtl) 

COMMON / LWPKV1 / IJVEC(500) 

DATA EPS / 1.E-2S / 

C 

C REARRANGE ROWS AND COLUMNS OF SMALL DENSE MATRICES ( A) EY SMALL IVECS 
C BOTH FROM NUT A , AND ADD MATRIX ELEMENTS AT LIKE LOCATIONS TO FORM 
C LARGE SPA* .> MATRIX (Z) ON NUTZ. NUT2 IS INITIATED IN THIS SUBROUTINE 
C TH.S, ANY r .<EVICUS DATA ON NUTZ IS DESTROYED. 

C MATRIX (A) ELEMENTS WHOSE ABSOLUTE VALUE IS LESS THAN EPS ARE NOT 
C PL/CED INTO MATRIX (Z). 

C A IS ASSUMED SYMMETRIC, ONLY THE LOWER HALF IS USEO. 

C ONLY THE LOWER HALF E Z (ALSO SYMMETRIC) IS FORMED. 

C CALLS FORMA SUBROUTINES YIN ,YINI ,YOUT ,YOUTI ,YPART ,YZERO , 

ZZBOMB. 

DEVELC'PFD FY R A PHILIPPUS. JANUARY 1973. 

LAST REVISION BY PL WOHLEN FOP. NASA. MAY 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH SMALL MATRIX DATA IS 
STORED. DATA CONSISTS OF SIZE, MATRIX ELEMENTS, IVEC. 

NUT A IS READ IN THIS SUBROUTINE WITH A READ STATEMENT. THUS, 
IT MUST HAVE BEEN GENERATED WITH A WRITE STATEMENT. NUTA 
CANNOT BE USED WITH ANY OTHER SPARSE SUBROUTINE. 

NUTZ = LOGICAL NUMBFR OF UTILITY TAPE ON WHICH MATRIX Z IS STORED, 
f NRZ = NUMPEP OF ROWS AND COLUMNS OF SPARSE MATRIX Z. 

i W = MATRIX WORK SPACE. 

C KW = ROW DIMENSION OF W IN CALLING PROGRAM. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPF. 

C NUT2 = LOGICAL NUMBER OF U T 1 LITY TAPE. 

C NUT3 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NERROR EXPLANATION 

C 1 = ROW OP COLUMN LOCATION OUTSIDE MATRIX Z. 

C 2 * NUTA CONTAINED ONLY NULL MATRICES. 

C 

CALL YZFRO (NUTZ , NRZ , NRZ ) 

REWIND NUTA 

READ (NUTA) NAME 

IF (NAMF.EC.6H ) RETURN 

REWIND NUT A 

REWIND NUT! 

REWIND NUTZ 
KVC4=KV/4 

C SET (2) MATRIX HEADER. 

MHFAD(l) - NRZ 
MHEAD ( ? ) = NRZ 
MHE AD ( 3 ) = 0 
MHEADCM = 0 
MHEAD(S) = 6HORDER 
MHEAD(c) = KV 


L ’T' 



YRVAD2 — 2/ 7 


MHEAD(7) = 5HLCWER 
ISHAP = 5HL0WER 
'C 

LZ = 0 
NSEG = C 
10 IRET = 1 

READ (NUTA) NAME,IE,N,(IP,I=1,7>,((W(I.J),I=1,N),J=I,N>, 

* (IJVEC(I) ,1=1,N) 

IF (NAME. EC. 6H ) GO TO 72 

NERRCR=1 

DO 15 1 = 1, N 

IF (IABS(IJVECd) J.GT.NRZ) GO TO 999 
15 CONTINUE 
IA = 0 

20 IA = TA+1 

IF (IA.GT.N) GO TO 10 
JA = 0 

30 JA = JA+1 

IF (JA.C-T.N) GO TO 20 
IF (IA.LT.JA) GO TO 20 
IF (ABS(W(IA,JA) J.LT.EPS) GO TO 30 
IF (IA.EO.JA) GO TO 35 

IF (IABS(IJVEC(IA)).NE.IABS(IJVEC< JAl)J GO TO 35 
W(IA,JA) = 2.*W(IA,JA) 

35 LZ = LZ*1 
TRET = 2 

IF (LZ.GT.KVCA) GO TO 72 
IF ( I JVEC ( I A ) ) 40,50,60 
40 IF ( I JVEC ( J A ) ) 45,50,55 
45 LV(LZ) = — I OCCOO*IJVEC ( IA I— I JVEC (JA I 
V ( L Z ) = W( TA,JA ) 

GO TO 30 
50 LZ = L7.-1 
GO TO 3C 

55 LV(LZ> = — TOCCOO*IJV£C (IA)+IJVEC( JA ) 

V(LZ) = — W ( I A , J A ) 

GO TO 30 

60 IF (I JVEC ( JA) I 65,50,70 

65 LV(LZ) = 1 CCOCC* IJVEC t IA)— 1 JVEC( JA ) 

V(LZ) = — W ( I A » J A ) 

GO TO 30 

70 LV(LZ) = ! OOOCC*IJVEC ( IA ) +1 JVEC ( JA ) 

V( LZ ) = W(IA,JA) 

GO TO SC- 

72 IF (LZ.GT.KVOA) L2=LZ-1 
IF (LZ.Cw.O) GO TO 225 
DO 65 K = 1 , .?. 

I = LV(K)/i OOCC 
J = LV(K) -1 00000*1 
IF (I.LT.J) tV(K)=100000*J+I 
85 CONTINUE 
r 

USE YLOPO LOGIC PUT ADD VALUES AT LIKE LOCATIONS. 

C SINGLETON METHOD 
M = 1 
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LAEM1 = L2-1 
1 = 1 

J = LZ 

105 IF (I.GE.J) GO TO 170 
110 K=1 

I J= ( J+I )/2 
IT=LVf IJ) 

IFdV(I).LE.IT) GOTO 120 
LV(IJ)=LVtI) 

LV(I)=IT 
IT=LV ( IJ) 

TG=V(IJ) 

V(1J)=V(I) 

V 1 1 )=TG 
120 L=J 

IFtlV( JJ.CE.IT) GO TO 140 
LV(IJ)'LVCJ) 

L V( J ) =IT 
IT=LV*IJ) 

TG- VIIJ ) 

v<ij>=vm 

V( J)=TC 

IF(LVm.LF.lT) GO TO 140 
LVlIJ)=TV<n 
L V( I ) = IT 

n^Lvtu) 

TC=V(IJ) 

v<Tj)=vm 

V!I )=TG 
GO TO 140 

130 LV(L)=LV(K) 

LV(K)-ITT 

TG=V(L) 

V(L)=V(K) 

VtK)=TG 
140 L — L — 1 

IF(LVIL).GT.IT) GOTO 140 
ITT=LV( L) 

150 K=K+1 

IF(LVfK I.UT.IT) GO TO 150 
IF(K.LE.L) GO TO 130 
IF(L-I.'-E.J-K) GO TO 16G 
I L ( M ) - I 
IU(M)=L 
I =K 
M=M+1 
GO TO 1E0 
160 1L (M)=K 
IU(M)=J 
J=L 
M=M + ] 

GO TO 1 80 
170 M=M-1 

IF(M.EC.O) GO TO 21C 
I=1L(M) 
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C 


c 


J=IU(M) 

18C IF(J-I.GE .11 ) GD TD 1 IC 
IF U.EC-.l) GC TC 105 
1 = 1-1 
19C 1=1+1 

IF(I.EC.J) GC TC 170 
IT=LV ( I+l ) 

IF ( LV ( I ) . LE. IT ) GC TO 19C 
TG=V (I + l ) 

K= I 

200 LV (K+l )=LV(K) 

V ( K+l ) =V ( K ) 

K=K— 1 

IF(IT.LT.LV(K)> GOTO 200 
LV ( K+l ) = IT 
V ( K+l ) =TG 
GC TC 190 


210 LZE = 1 
LAW = 1 

IF (LZ.EQ.l) GC TO 22 2 
DC 22C 1=2, LZ 

IF ( LV (LAW ) .EC. LV( I ) ) V(t AW1=V(LAWJ+V( II 
VCLZE )=V(LAW) 

LV(LZE)-LV(LAW) 

IE (LV(LAW).EQ.LV(l)) GC TO 220 

LZE=LZE+1 

LAW=I 

IF (I.LT.LZ) GO TO 220 
V ( LZ E ) = V( I ) 

LV(LZEf=LV( I > 

220 CONTINUF 
222 NS EG = NSEG+I 
MPHEAD(l) = LZE 
MPHEAD ( 2 ) = L V ( 1 ) 

MPHEAP ( 3 ) = LV(LZF) 

CALL YCL’T I (NUT1 , MPHEAD, 1 ,10) 

CALL YOUTI (NUT I ,LV ,1 , L ZE ) 

CALL YCUT (NUT1 , V,1,LZE> 

LZ = 0 

GO TO (225,35), 1RET 


225 IF (NSEG.GT . 1 ) GC TC 228 


IF (NSEG.tT.l) GO TD 999 
REWIND NUTZ 
MHEAD ( 3 ) = 1 
MHE AD (4) = LZE 

CALL YCUT I (NUTZ »MHE AD ,1,10) 
MPHEAD ( 1 ) = LZ F 
MPHFAD ( 2 ) = LV ( 1 ) 

MPHEADO) = LV(LZE) 

CALL YOUTI (NUTZ , MPHEAD, 1, 10) 
CALL YOUTI (NUTZ, LV,1, LZE) 
CALL YCUT (NUTZ, V,1,LZE) 


NERR0R=2 
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GO TO 410 

C NOW THERE ARE NSEG ORDERED GROUPS WRITTEN ON NUT1. 
22 8 NT1 = NUT 1 
NT2=NUT2 
NPFC = NSEG 
NSEG = 0 
NNZ2 = C 
REWIND NUT3 

MESHING OPERATION 
23C REWIND NTI 
RFWIND NT2 

IF (N^'EC.EC.D GO TO 305 
CALL T INI (NTI, MPHEAD, 1,10) 

NN2P1 = MPHEADd ) 

CALL YINI (NTI , LV, 1 , NNZP 1 ) 

CALL YIN (NTI ,V,1,NNZP1) 

IF (NPEC.EQ.l) GO TO 305 
NPFC2=0 
L2E = NNZP1 

1 = 1 

235 I = 1+1 

IF (I.GT.NREC) GO TO 300 
CALL YINI (NTI ,M2HEAD, 1,10) 

NN2P2 = M?HEAD(1 ) 

LP2S = NNZP1+1 
LP2E=LP?S— 5 +NNZP2 
CALL YINI (NTI »IV,LP2S,LP2E) 

CALL YIN (N T 1 ,V,LP ?S, LP2E ) 

IF ( LV ( LP2S ) .GT. LV(NNZPD) GC TO 291 

MESH TWO PARTITIONS 
11=1 

I? = LP2S 
IW=2*KV04 
IZ=0 

250 IW=IW+1 

IF (LV(I1)-LV(I2) ) 265,265,255 
255 V ( IW )=V (12) 

L V( IW )=LV ( 12 ) 

12 = 12+1 

IF (I2.GT.LP2E) GO TO 275 
GO TO 250 
265 V(IW)=V(I!> 

LV(IW)=LV(11) 

11 = 11+1 

IF (I1.GT.LP2S-1 ) GC TO 285 
GO TO 250 

275 NELTM = LP2S-I1 
K=LP2E 
L = IP2S-! 


DO 280 J=1,NFLTM 



YRVAD2 — fc/ 7 


V(K )=V( L) 

LV ( K )=LV( L ) 

K=K— 1 
28C L-L-l 
C 

285 IF C IW.FQ.2*KV04) GO TO 29! 

J1=2*KV04+1 

c 

DO 290 J=J1,IW 
IZ=IZ+1 
V( JZ)=V(J ) 

290 LV ( IZ ) = LV ( J ) 

291 LZE = I 
LAW = 1 

DO 293 J=2»LP2E 

IF (LV(LAW) .FG.LVU)) V(LAW)=V(LAW) +V(J ) 
LV(LZF) = LV (LAW ) 

V(LZF) = V(LAW J 
IF (LV(LAW) .EQ.LV(J) ) GO TO 293 
LZE = LZE +1 
LAW - J 

IF (J.LT.LP2E) GO TO 293 
LV(LZF) = LV ( J ) 

V ( LZE ) = V(J) 

293 CONTINUE 

IF (LZE .LT.KVO** .AND. NREC2.EG.C) NNZP1=LZE 
IF (LZE.tT.KV04 . AND. NREC2.EQ.0 1 GO TO 235 
IF (LZE.E0.NNZP1) GO TO 235 
C 

NR EC 2 = NREC2+1 

M2HE AD ( I I = LZE-NNZP1 

M2HEAD(?1 = LVCNNZP1 + 1) 

M2HEAD ( 3 ) = LV(LZF) 

CALL YOU! I (NT2,M2HFAD,1,1C) 

CALL YCUTI (NT2 ,LV,NNZP1+1 ,LZE) 

CALL YCUT (NT2, V,NNZP1+1,LZE) 

GO TO 235 
C 

C ALL NRFC PARTITIONS HAVE BEEN READ FROM NT1 
300 MPHFAD(l) = NNZP1 
MPHEAD ( 2 ) = LVI]} 

MPHFAD ( 3 ) = LV(NNZPl) 

CALL YCUTI (NUT 3 »MPHE AD »I ? IC ? 

CALL YCUTI (NU7 3 t LV f 1 ,NNZP1> 

CALL YOUT (NUT3, V,1,NNZP1> 

NSEG = NSFG+J 
NNZZ = NNZZ+NNZP1 
NRFC=VPFC2 
NTS=NT) 

NT1=NT2 
NT2=NTS 
GO TO 230 

305 IF (NPFC.EQ.C) GO TO 400 

CALL YOUTI (NUT 3 »MPh E AD » 1 1 10 ) 
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CALL YOUTI (NUT3.LV.1.NNZP1 ) 
CALL YOL'T INUT3, V.1.NN2P1) 
NSEC = NSFG+l 
NNZZ = NNZZ+NN2P1 
400 REWIND NUT3 
REWIND NUTZ 
MHE AD ( 3 ) = NSEC 
MHEAD ( A ) = NNZZ 
CALL YOUTI (NUTZ .MHEAD.l.lO) 

DO 4C5 1= l.NSEG 

CALL YTNI ( NUT 3 ,MH E AD , 1 , 1 0 ) 

CALL YOUTI (NUTZ .MHEAD,! .10) 

K = MHEAD ( I ) 

CALL YINI (NUT3.LV, l.K) 

CALL YOUTI (NUTZ ,LV,! ,K) 

CALL YIN (NUT3, V.l.K) 

405 CALL YOUT (NUTZ, V.l.K) 

410 CALL YN0ZER (NUTZ .V.LV.KV.MJTl ) 
RETURN 

999 CALL ZZEOMB (6HYR VAD2 .NERROR ) 
END 
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SUBROUTINE YRVAD3 ( NUTA ,NUTZ ,NRZ t NCZ , W, KW , V,LV ,KV ,NUT 1 ,NUT2 t NUT3 ) 
DIMENSION VC11 ,L V (1 ) , IU < 16 ) , IL < 16 ) , MHEAP ( 1 D) ,MPHE AD( 10 ) ,M2HEAD< 1C ] 
DIMFNf TON W ( KW , 1 ) 

COMMON / LWPKV1 / I VECI25C), JVEC(25G) 

DATA EPS / l.E-25 / 

C 

C REARRANGE ROWS AND COLUMNS OF SMALL DENSE MATRICES (A) BY SMALL 1VEC! 
C AND JVECS, ALL FPL'M NUTA, AND ADD MATP IX ELEMENTS AT LIKE LOCATIONS 
C TO FORM LARGE SPARSE MATRIX ( 2 ) ON NUTZ. NUT 2 IS INITIATED IN THIS 
C SUBROUTINE. THUS, ANY PREVIOUS DATA CN NUTZ ARE DESTROYED. 

C MATRIX (A) ELEMENTS WHOSE ABSOLUTE VALUE IS LESS THAN EPS ARE NOT 
C PLACED INTO MATRIX € Z) . (A) IS SQUARE. 

C CALLS FORMA SUBROUTINES YIN ,YIN1 ,YOUT ,YOUTI , YPARI ,2ZB0M6. 

C DEVELOPED BY R A PHILIPPUS. FEBRUARY I97S. 

C LAST REVISION BY WA bENFIELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS ( ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH SMALL MATRIX DATA IS 
C STOPFD. DATA CONSISTS OF SIZF, MATRIX ELEMENTS, IVEC. 

C NUTA IS READ IN THIS SUBROUTINE WITH A READ STATEMENT. THUS, 

C IT MUST HAVE BEEN GENERATED WITH A WRITE STATEMENT. NUTA 

C CANNOT BE USED WITH ANY OTHER SPARSE SUBROUTINE. 

C NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX 2 IS STORED. 

C NFZ = NUMBER OF ROWS CF SPARSE MATRIX Z. 

C NCZ = NUMBER OF COLUMNS OF SPARSF MATRIX Z. 

C W = MATRIX WORK SPACE. 

j* KW = ROW DIMENSION OF W IN CALLING PROGRAM. 

I V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT3 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NEPPOF EXPLANATION 

C 1 = LOCAL SIZE EXCEEDS FINAL SIZE. 

C 2 = NUMBFR OF SEGMENTS LESS THAN ONE. 

C 

REWIND NUTA 
REWIND NUT 1 
RFWIND NUTZ 
KV04=KV/4 
MHE AD( 1 ) - NPZ 
MHE AD ( 2 ) = Nr,z 
MHFAD(S) = ■ ’ORDER 

MHFAD ( 6 ) = C 
MHEAD( 7 ) = 6HWH0LE 
ISHAP - fcHWHOLE 
C 

L Z = 0 
NSEG = 0 
10 IRET = I 

RFAD (NUTA) NAME , IP ,N , ( IB , 1=1 ,7) , ( ( W ( I , J ) , 1=1 ,N) , J=1 ,N ) , 

* (I VEC(I) ,I=1,N), ( JVEC ( I ) , 1=1 *N) 

IF (NAME . EQ.6H ) GO TO 72 



o o 


YRVAD3 — ?/ 7 


NERR0R=1 

DO 15 1-1 *N 

IF (IAFS(I VEC(I ) l.GT.NRZ) GO TO 999 
IF (IAP5(J VEC(I ) J.GT.NCZ) GO TO 999 
15 CONTINUE 
IA = 0 

20 IA = I A+l 

IF (1A.GT.N) GO TO 10 
JA = 0 

30 JA = JA+1 

IF (JA.GT.N) GO TO 20 
IF (ABE (W (IA»JA) ) .LT.EPS) GO TO 30 
35 LZ = LZ+1 
IRET = 2 

IF (LZ.C-T.KV04) CO TO 72 
IF (I VEC(IM) 40,50,60 
40 IF ( JVEC(JA)) 45,50,55 
45 LV(LZ) = -100000*1 VEC ( I A )— JVECiJA) 

V ( LZ ) = W ( IA* J A ) 

GO TO 3C 
50 LZ = LZ-I 
GO TO 3C 

55 LV(LZ) = — 1 OGOOO*I VEC(IA)+ JVEC(JA) 

V 1 LZ ) = — W ( I A , J A ) 

GO TO 30 

60 IF ( JVEC(JA)) t>5 ,50, 70 

65 LV(LZ) = 100000*1 VEC(IA)- JVFC(JA) 

V ( LZ ) = -W(IA,JA) 

GO TO 30 

70 LV(LZ) = 1C0C00*I VEC(IA)+ JVEC(JA) 

V( LZ ) a W ( I A , J A ) 

GO TO 3C 

72 IF (IZ.GT.KV04) LZ=LZ-1 
IF (LZ.EQ.O) GO TO 225 

SINGLETON METHOD 
M = 1 

LAFM1 = LZ— 1 
1 = 1 

J = LZ 

105 IF (I.CE.J) GO TO 170 
110 K= I 

I J= ( J+I )/2 
IT-LV(IJ) 

IF ( LV( D.LE.IT) GO TO 120 
LV ( I J ) =LV (I ) 

L V( I ) =IT 
ITdV( IJ) 

T G=V ( I J ) 

Vd J)=V(I ) 

V (1 * =TG 
120 L = J 

IFILV(J).CE.IT) GO TO 140 
LV( IJ)=LV( J) 

LV( J ) -IT 
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1T=LV ( I J ) 

TG=VCIJ) 

V(IJ)=V(J) 

V ( J )=TG 

IF ( LV ( I ) • LF. IT) GOTO 140 
LV( T J ) = LV (I ) 

LV(1I=IT 
IT=LV< I J) 

TG=VUJ) 

V 1 1 J )=V C I I 
V( 1 J =TG 
CO TP 140 

130 mn=LV(K) 

LV(K)=ITT 

TG=V(L) 

V(L)=VCK) 

V (K )=TG 
140 1=1-1 

IFCLVCLI.GT.ITI 60 TO 140 
ITT=IV<1) 

150 K=K + T 

IF(LVIK).LT.IT) GO TO 150 
IF(K.LE.L) f TP 130 
IFtL-I.LE.J-K) GO TO 160 
1 l ( M )=I 
IU ( M } 

1= K 
M-- M+J 
GO TP 180 
160 IL(M)=K 
IU(M)=J 
J=L 
M=M+1 
GO TO 180 
170 M=M-1 

IFIM.FC.O) GO TO 210 
I-IL I M) 

J=IU(M) 

180 IFCJ-I.GE.Jl) GO TO 110 
IF (l.EO.n GO TO 105 
1 = 1-1 
190 1=1+1 

IF(J.FQ.J) GO TO 170 
1T=LV ( 1 + 1 ) 

IF(LVU I.LE.IT1 GO TO 190 
TG=V <1 + 1 ) 

K = 1 

200 L V ( K+l ) =L VI K ) 

VlK + 1 )=V( K ) 

K = K— 1 

IFCIT.LT. LVIKH GO TO 200 
L V I K+l ) = 1 7 

V ( K+l ) = TG 

GO TO 190 


C 
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210 LZE = 1 
LAW - 1 

IF (L2.F0.1 ) GO 70 222 
DC 220 1=2, L2 

IF(LV(LAW) .FQ. LV(IU V( LAW ) s V( LAW ) +V( I ) 
V ( LZF ) =V ( LAW ) 

LVILZF )=LV(LAW ) 

IF (LV(LAW) .EQ.LV(I) ) GO TO 220 
LZE=L7.E + 1 
L AW= I 

IF (I.LT.LZ) GO TO 220 
V(IZF:)=V(I) 

LV( LZF ) =LV( I ) 

220 CONTINUE 
222 NSEG = NSFG+1 
MPHEAD(l) = LZF 
MPHEAD ( 2 ) = L V ( 1 ) 

MPHEAD ( 3 ) = LV ( L 7. E) 

CALL YOUTT ( NUT1 ,MPHE AD , * , 10 ) 

CALL YOUTI (NUT1 ,LV,1 ,LZt > 

CALL YOUT (NUT1 , V,1,LZE) 

LZ = C 

GO TO (225,35) , IRET 

225 IF (NSEG.GT.l) GO TO 228 

IF (NSEP.LT.l) GO TO 999 
REWIND NUTZ 
MHEAD(3 ) = 1 
MHEAD(4) = LZE 

CALL YOUTI (NUTZ ,MHE AD , 1 , 10 ) 

MPHEAD(l) = LZF 
MPHEAD12) = LV ( I ) 

MPHEAD(3> = LV(LZF) 

CALL YOUTI (NUTZ , MPHEAD, 1,10) 

CALL YOUTI (NUTZ ,LV,1, LZF) 

CALL YOUT (NUTZ, V,1,LZE) 

GO TO 410 


NERRCR=2 


NOW THERE ARE NSEG ORDERED GROUPS WRITTEN ON NUT1. 
228 NT1 = NUT1 
NT2=NUT2 
NR EC = NSEG 
NSFG = C 
NNZZ = 0 
REWIND NUT3 


MESHING OPERATION 
230 REWIND NT1 
REWIND NT 2 

IF (NRFC.EQ.C) GO TO 305 
CALL Y1NI (Ni l , MPHE AD 't » 1C ) 
NNZP1 = MPHEADU ) 

CALL YINI (NT1 ,LV,1,NNZ" 1) 

CALL YIN (NT) ,V,1 ,NNZP1 ) 
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IF (NRFC.EG.l ) GO TO 305 
NR E C2=0 
LZE = NNZP1 
C 

I = I 

235 I = 1 + 1 

IF (I.GT.NREC) GO TO 300 
CALL YINI (NT1 ,M2HEAD,1,1G) 

NNZP2 = M2HE AD ( 1 ) 

LP2S = NNZP1+1 
LP2E=LP2S-1 +NNZP2 
CALL YINI (NT1 , LV , L P2S , LP2E ) 

CALL YIN (NT! »V, LP2S»LP2F) 

IF ( LV( LP2S ) .6T. LV(NNZPD) GO TO 291 

MESH TWO PARTITIONS 

II = 1 

12 = LP2S 
IW=2*KV0A 
IZ=0 

250 IW=1W+1 

IF (LV(I1 >-LV(I2)» 265,265,255 
255 V ( IW )=V ( 12 ) 

LV ( IW) =LV ( 12 } 

12=12+1 

IF (I2.GT.LP2E) GO TO 275 
GO TO 250 
265 V(IWi=V(II) 

LV ( IW ) =LV ( 1 1 ) 

11=11+1 

IF U1.GT.IP2S-1) GO TO 285 
GO TO 250 

275 NFLTM = LP2S-II 
K=LP2F 
L = LP2 $~1 

DO 280 J=1 ,NELTM 
V ( K ) =V ( L ) 

LV( K )=LV( L ) 

K = K-1 
280 L=L-1 

285 IF (IW.F0.2*KVOA) GO TO 291 
J1=2*KV04+1 

DO 290 J=J1,IW 
IZ=1Z+1 
V(IZ)=V(J) 

290 LV(1Z)=LV(J) 

291 LZF = 1 
LAW = 1 

DO 293 J=2,LP2E 

IF (LV(LAW) .FQ.LVU) ) V ( LAW ) =V ( LAW ) +V ( J 1 
LV(LZE) = LV ( LAW ) 

V(LZE) = V ( LAW ) 
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IF (LV(LAW).EQ.LV(J)) CD ID 293 
LZF = LZE + 1 
LAW = J 

IF (J.LT.LP2E) GO TO 293 
LV(LZE) - LV ( J ) 

V(LZF) = V(J) 

293 CONTINUE 

IF <LZE.LT.KV04 .AND. NREC2.EC.0) NNZPl^LZF 
IF (LZF.LT.KV04 . AND. NREC2.EG .0 ) GO TO 235 
IF (LZF.FQ.NNZP1 J GO TO 235 

NRFC2 = NRFC2+1 
M2HFAD(1) =• LZF-NNZP1 
M2HE AD ( 2 ) = LV(NN7P1+1) 

M2HE AD ( 3 ) = LV(LZF) 

CALL YOUTI (NT2,M2HEAD,1,10) 

CALL YOUTI (N"»2,LV,NNZP1 + 1,LZE) 

CALL YOUT ( NT2 » V,NNZP1+1,LZE) 

GO T^ 235 

ALL MRFC PARTITIONS HAVE BEEN READ FROM NT1 
300 MPHEAD(l) = NNZP1 
MPHEAD (?) = LVU) 

MPHEAD ( 3) = LV(NNZPI) 

CALL YOUTI (NUT3, MPHEAD, 1,10) 

CALL YOUTI (NUT3 ,LV, 1 ,NNZP! ) 

CALL YOUT (NUT3, V,1,NNZP1) 

NS EG = NSFG+1 

NNZZ = NNZ Z+NNZP1 

NREC=NREC2 

NT $=NT 1 

NT1=NT2 

NT2-NTP 

GO TO ... 0 

3C5 IF (NKEC.FQ.C) GO TO 400 

CALL Y' JTI ( NUT 3 ,MPH F AD , 1 » 1C ) 

CALL YOUTI (NUT3 ,LV, 1 ,NNZP1 ) 

CALL YOUT <NUT3, V»1,NNZP1) 

NS EG = NSFG+1 
NNZZ = NNZZ+NNZPl 
400 REWIND NUT 3 
REWIND NUTZ 
MHEADC3I = NSEG 
MHE AD ( 4 ) = NN? 7 
CALL YOUTI (NUTZ »MHE AD, I » 10 ) 

DC 405 T=1,NSFG 

CALL Y INI (NUT3,MHEAD,1,1G) 

CALL YOUTI (NUTZ , MHE AD, 1,10) 

K - MHFAD(l) 

CALL Y1N1 (NUT 3 ,LV, 1 ,K ) 

CALL YOUTI (NUT 2 ,LV, 1 »K ) 

CALL YIN (NUTT, V,1,K) 

405 CALL YOUT (NUTZ, V,!,K) 

410 CALL YNOZFR (NUTZ ,V,L V ,KV,NUT1 ) 
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RETURN 

999 CALL Z 2 BOMB (6HYRVAD3 »NERROR ) 
END 
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SUBROUTINE YRVISI C A, JVEC ,NU72 ,NP A2 *NCA ,NC Z,V, LV, KV.KA ) 

DIMENSION V(l), LV(1), JVEC(l), MHEAD(IO), A(KA,l) 

DATA EPS / l.E-25 / 

REARRANGE COLUMNS OF DENSE MATRIX A TO FORM SPARSE MATRIX Z. 

ROWS OF A ARE NOT REARRANGED. 

CALLS FORMA SUBROUTINES XLORD ,YOUT ,YOUTI ,ZZBCMB. 

DEVELOPED BY R A PHILIPPUS. JANUARY 1973. 

LAST REVISION BY k'A BENFIELD. MARCH 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

A = MATRIX A. 

JVEC = VECTOR (SIZE = NCA ) . 

JVECIJ) = COLUMN POSITION OF A(COL J) IN Z- 
IP JVFC(J) IS PLUS , 2 = A (COL J). 

IF JVEC(J) IS MINUS, Z = -A(COL J). 

IF JVEC(J) IS ZERO , A( COL J) IS OMITTED IN Z. 

NUTZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX Z IS STORED. 
NRAZ = NUMBER OF ROWS OF A,Z. 

NCA = NUMBER OF COLUMNS OF A. 

NCZ - NUMBER OF COLUMNS OF Z. 

V = VFCTCP WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SI2E OF V,LV IN CALLING PROGRAM. 

KA = ROW DIMENSION OF MATRIX A IN CALLING PROGRAM. 


NERRCR EXPLANATION 

1 = COLUMN LOCATION OUTSIDE MATRIX Z. 

2 = DIMENSION SIZE EXCEEDED (KV). 


REWIND NUTZ 
C 

NERRCR=1 

DO 2 1-1 ,NCA 

IF ( I ABS ( JVEC( I ) ) .GT.NCZ ) GO TO 999 
2 CONTINUE 
C 

J = 0 

NERROR=2 

DO SO I A = 1 »NR AZ 
IA1 = 1 CCCCC*I A 
DO 50 JA=!,NCA 

IF (APSIA(IA,JA) ) .LT.EPS) GO TO 50 
J = J+l 



IF (J. 

GT.KV) 

GC TC 999 


IF (JVFC(JA) 

) 40,25,45 

25 

J - J- 

1 



GO Tn 

50 


AO 

V(J) = 

—A ( I A 

, JA ) 


LV(J) 

= I Al- 

JVEC ( JA ) 


GO TO 

50 


45 

L V ( J i 

= IAI + 

JVEC ( JA) 


V(J) = 

A ( 1A » 

JA) 

50 

CONTINUE 



C 
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MHE ADf 1 ) = NRAZ 

MHFAD (2 ) = NCZ 

MHE AD (3 ) = 1 

MHEAD(A) = J 

MHE AD ( 5 ) = 5HCRDER 

MHE AD ( 6 ) = KV 

MHEAD ( 7 ) = 5HWH0LE 

IF (J.GT.KV/4) MHEAD ( 6 ) - 0 

CALL XL CRD (V,LV,1,J) 

CALL YCUTI (NUTZ * MHEAD, 1,10 
MHE AD (I) = J 
MHEAD12) =■ LV Cl) 

MHE ADC 3 ) = LVC J ) 

MHEAD ( 4 ) = 0 
MHEAD 15} = 0 
MHE AD (6 ) = 0 
MHEAD 17) =■ 0 

CALL YCUTI (NUTZ,MHEAD t l,!0) 
CALL YCUTI (NUTZ,LV»I tJ) 

CALL YCUT (NUTZ »V* I » J ) 

C 

RETURN 

999 CALL ZZBCMB (6HYRVIS1 *NERR0R ) 
END 
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SUBROUTINE YRVTOD (NUT A , IVEC, JVEC ,2 , NRZ ,NCZ ,V, LV,K V,KR2 ) 
DIMENSION IVEC (I } ,JVEC ( 1 ) , Z ( KRZ, 1 ) , V( 1 ) ,LV (1) ,MH ( 10 ) 

REARRANGE AND ADD ROWS AND COLUMNS OF SPARSE MATRIX A 
TO DENSE MATRIX 2. 

PE SURE TO DEFINE 2 BEFORE CALLING THIS SUBROUTINE. 

FOR EXAMPLE CALL SUBROUTINE ZERO TO SET MATRIX 2 TO ZERO. 
CALLS FORMA SUBROUTINES YINJ »YIN AND ZZBOMB. 

CODED BY JOHN ADMIRE *NASA* DECEMBER 1®73 . 


C 

C 


SUBROUTINE ARGUMENTS 

NUTA - INPUT LOGICAL NUMBER OF UTILITY TAPE ON WHICH 
SPARSE MATRIX A IS STORED. 

IVEC - INPUT I VEC ( 1 ) = ROW POSITION OE A(ROW I) IN 2. 

IF IVEC ( I ) IS PLUS ,2 =Z (ROW IVEC (I ) )+A(POW 1). 

IF IVEC ( I ) IS MI NUS ,Z =Z (ROW IVEC (1 ) )-A (ROW I). 

IE IVEC (I) IS ZERO , ROW I OE A COES NOT APPEAR IN 2. 
JVEC - INPUT JVEC(J)= COL POSITION OF A(COL J) IN 2. 

IF JVEC(J) IS PLUS *2=2 (COL JVEC(J) )+A(COL J). 

IF JVEC(J) IS Ml NUS » 2 =Z (COL JVEC (J 1 >— A ( COL J). 

IF JVEC(J) IS ZERO ,CCL J OF A COES NOT APPEAR IN 2. 
2 - INPUT/OUTPUT MATRIX TO WHICH ELEMENTS OF A ARE ADDED OR 

SUBSTRACTED. 

NR2 - INPUT NUMBER OE ROWS IN Z. 

NCZ - INPUT NUMBER OF COLS IN 2. 

V - VECTOR WORK SPACE. 

LV - VECTOR WORK SPACE. 

KV - INPUT DIMENSION OF V AND LV IN CALLING PROGRAM. 

KRZ - INPUT DIMENSION OF NUMBER OF ROWS OF 2 IN CALLING PROGRAM. 


NERR0R=1 

I F ( NR2 .GT. KRZ) GO TO 999 
REWIND NUTA 

CALL YINI (NUTA, MH, 1,10) 
NRA=MH( 1 ) 

NCA=MHf 2 ) 

NPART=MH(3) 

IF(NPART .EO. C) return 


NB =0 


DO 1C 1 = 1 ,NR A 

IF( 1AES (IVEC (I ) ) .GT. NB ) NB = I ABS ( I VEC ( I ) ) 
10 CONTINUE 
NERR0R=2 

IF (NB .GT. NRZ ) GO TO 999 
NB = 0 

DC 20 J = 1 , NC A 

IF ( I APS (JVEC ( J ) ) .GT. NB) NB = T ABS ( JVEC ( J ? ) 
20 CONTINUE 


NERR0R=3 

I F (NB .GT. NC2) GC TO 999 
DO 110 K=1 ,NPAPT 
CALL YINI (NUTA, MH, 1,10) 
NNZP=MH(1 ) 

CALL YINI (NUTA, LV,1,NNZP) 
CALL YIN (NUTA,V,1,NNZP) 
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DO 100 ll-I»NNZP 
I=LV(LL)/1CCC0C 
IFIIVFCCI ) ) 3 C » 1 C C *40 
III. II = IABS(IVE.C(I)) 

NS = -1 
GO TC 50 
■o ii=ivec(I) 

NS=1 

£0 J=LV f LL )— 1 GOCCC* I 
IF! JVEC( J))60,1CC t 7C 
60 JJ=IAr>S(JVEC(JJi 
NS=-NS 
GO TO 80 
T ? JJ= JVEC ( J ) 
t ) IF (NS .LT. C) GO TO 9C 
Z(II*JJ)=Z(IItJJ )+V(LL) 

GO TO 100 

90 Z( II,JJ)=Z(II,JJ)-V(LL> 

100 CONTINUE 

110 CONTINUE 
RETURN 

999 CALL ZZ.<GM6 (6HYR VTOD, NERROR } 
END 
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SUBROUTINE YSRED2 (NUTA ,NUTR ,NUTT,NR , IFT, V ,LV ,KV ,NUI 1 *NUT2,NUT3 * 

* NUT4I 

DIMENSION V< 1),LV(1) , MHEAD(IO) ,MPHEAD(10) , IhC 10) 

DATA EPS/1. E-2C/ 

REDUCE SPARSE STIFFNESS MATRIX (A) TO FORM SPARSE REDUCED STIFFNESS 
MATRIX <P) AND ION OPTION) REDUCTION TRANS FORMATION MATRIX (T). 
COORDINATES TO BF RETAINED APE NUMBERED LAST. 

IF THE WHOLE MATRIX (A) IS INPUT, ONLY THE LOWER HALF IS USED. 

RAND WIDTH (DIAGONAL UP TO TCP NON— ZERO) MUST BE LESS THAN OR EQUAL 
TO I KV— N ) / 2 , WHFRE N IS MATRIX SIZE (SQUARE). 

CALLS FORMA SUBROUTINES YIN ,YINI ,YLORD ,YOuT ,YOUTI ,YPAR1 , 

YTPANS, ZZEOMB . 

DEVELOPED BY P L WOHLEN AND W A BENEIFLD. DECEMBER 1972*- 
LAST REVISION BY WA BENFIELD. MARCH 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH (A) IS STORED. 

NUTR = LOGICAL NUMEFR OF UTILITY TAPE ON WHICH ( R ) IS STORED. 

NUTT = LOGICAL NUMB FP OF UTI LT T Y TAPE CN WHICH (T) IS STORED. 

NR = NUMBER OF ROWS IN THE REDUCED STIFFNESS MATRIX. 

IFT = 0, TRANSFORMATION MATRIX (T) WTLL NOT PR CALCULATED. NUTT 
NEED NOT BE DEFINED IN CALLING PROGRAM. 

- 1, TRANSFORMATION MATRIX (T) WILL BE CALCULATED. 

V = VECTOR WORK SPACE. 

LV = VECTOR WOPK SPACE. 

, KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM, 

s. NUT1 = LOGICAL NUMBER 0* UTILITY TAPE. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT3 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE. 

C 

C NERROR E XPLANAT ION 

C 1 = BANDWIDTH LIMITATION EXCEEDFD (KV). 

C 2 - DIMENSION SIZE EXCEEDED (KV). 

C 3 = BANDWIDTH LIMITATION EXCEEDED (KV). 

C A = MATRIX IS SINGULAR. 

C 

C CONVERT (A) FROM SPARSE (NUTA) TO BAND (NUTI ) NOTATICN. 

KVP4=KV/A 
KV02 = KV/2 
KVC2R1 = KV02+1 
LAS=KV04+1 
REWIND NUTA 

CALL Y INI (NUT A ,MHE AD ,1,10 

NR A = MHEADU) 

NO = NR A - NR 
KVMN = KV-NRA 
KVMN02 = KVMN/? 

IASHAP - MHF AD (7 ) 

NUTS=NUTA S 

IE (IASUAE.EC.BHUPPER ) CALL YTR’ANS (NUTA, NUTR ,V,LV,KV,N'JT1 ,NUT2 ) 

IE ( IASHAP. EQ. 5HUPPER ) NUTS -NUTR 
CALL YLOPD (NUTS ,V,L V,KV* NUI 1 ,NUT? ) 

REWIND NUTS 


S 
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CALL YINI <NUTS,MHE AD, 1,10) 

REWIND NUT1 

ILV = KVC‘4 

JLV = KVOA+NRA 

IF ( JLV.LT.KV02) JLV=KV02 

JLVS = JLV 

KP = 1 

LAMAX = LA5-1+KVMN02 

LAE = LAS 

JS = 1 

NGROUP = 0 

LAS1 = KVC4 

00 5 I=LAS,KV 

LVI I ) = 0 

5 vm=c. 

NNZZ = C 

NPARTA = MHEADI3) 

NROWS - 1 

DO 20 1=1, NPARTA 

CALL YINI (NUTS,MPHEAD, 1,101 

NNZPA = MPHEADU) 

CALL YINI (NUTS ,LV, 1 , NNZPA) 

CALL YIN (NUTS,V,1 , NNZPA) 

00 ?C J=l, NNZPA 
IA=LV( J)/10000C 
JA=LV( J)-1 00000*1 A 
IF (1A.LT.JA) CO TO 2C 
IF (IA.FQ.KP) GO TO 15 
LAS 1 = LAE 
LAE = LAF+TA-JA+1 
NFLR = KP-JS+1 

NERRCR=1 

IF (NELR.GT.KVMN02) GO TO 999 
NNZZ = NNZZ+NELR 
KP = KP+1 
JS = JA 

NRCWS = NROWS+1 

ILV = I LV+1 

LV(ILV) = NFLR 

IF (LAF.LE. LAMAX) GO TO 15 

JLV = JLV+1 

NERRCR=2 

IF (JLV.GT.KV) GO TO 999 

NROWS = NROWS-1 

LV(JLV) = NROWS 

NROWS = 1 

LAE = LAE-IA+Ja-1 

NGROUP = NGROUP+1 

CALL YOUT (NUT1 ,V,LAS,LAt) 

DO 10 L.=LAS »LAE 

10 V ( L ) =0. 

L AS 1 = KVOA 

LAF = KVC4+1A-JA+1 

KP = IA 

15 LA = LASI+JA-JS+1 
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V ( L A ) =V ( J ) 

20 CONTINUE 

IF (LAS.GT.I AE) GO TO 65 
NGROUP = NGRCUP+ 1 
ILV = ILV+I 
LV(lLV) = KP-JS+1 

NERROR= 

IF IL V! 1LV) .GT.KVMNC2 ) GO TO 999 
NNZZ = NNZZ.+LV( 1 LV) 

•LV = JLV + 1 

IF tJLV.GT.KV) GC TO 999 

L V ( J LV ) = NROWS 

CALL YOUT (NUT1 ,V,LAS»LAE ) 

65 DO 30 1=1 ,NRA 
30 LVC i) = LVIKVOA + I ) 

DO 40 1=1, NGROUP 
40 LV ( KVC?«-I ) = LV ( JLVS + I ) 

C 

C REDUCTION. 

C D IN VII THRU N ) . A,U GROUP A START AT V(N+1). 

C A ,U GROUP E START AT VlN + l-MKV-NJ/2 ). 

C LV(I),1=1,N IS NUMBER OF ELEMENTS IN COLUMN I. 

C LV(KV/?+lG) IS NUMBER OF COLUMNS IN GROUP IG. 

N = NR A 
NG = NGROUP 
LSGA = N+l 

LSGB = LSGA + (KV-NJ/2 

REWIND NUT3 

JEGA = 0 

DO 1<?5 TO A= 1 ,NG 

REWIND NUT1 

REWIND NUT? 

NUTP = NUT! 

NUTQ = NUT? 

IF ( 2* C IGA/2 ) .FG. IGA) NUTP=NUT2 
IF (NUTP .FQ. NUT?) NUT0=NUT1 

C OPERATE ON GROUP A ONLY. 

NCGA = LV (KVC2+I GA) 

JSG A = JEGA+I 
JEGA = JSGA+NCGA-1 
LEGA = LSGA— 1 
DO 101 J= JSGA, JEGA 
101 LFGA = LEGA ♦ LV(J) 

CALL YIN (NUTP, V, LSGA, LEGA) 

LJJ - LSGA—1 

DO 14-0 J= JSGA, JEGA 

JM1 = J-l 

IFND = J - 1 

IF (J .GT. ND) I E ND= J 

ITPPJ = J-LV ( J ) + 1 

LITDPJ = LJJ+1 

LJJ = LITOPJ+LV( J )-l 

IF (J .FQ. JSGA .AND. J .LE. ND) GC TO 134 
IF ( I TOP J .FQ. J) GO TO 134 
ISTART = ITOPJ 


3/1C 


3 


1 



YSRED? — A/10 


LI J = LITOPJ-1 

IF (ITOPJ .GE. JSGA) GO TO 105 
I START = JSGA 
LI J = L 1 TOP J+ JSGA— ITOPJ— 1 
105 LITOPI - LSGA 

IF (ISTART .EO. JSGA) GO TO 110 
ISM1 = ISTART -1 
DO 107 I=JSGA» 1 SMI 
107 LITOPI = LITCP I+L V( I ) 

110 DO 128 I=ISTART,IEND 
KEND =1-1 

IF (I .GT. ND) KEND = ND 
L I J = LIJ+1 
S = V(LIJ) 

IM1 = 1-1 

ITOPI = I-LVU> + 1 

IF (ITC'PI .LT. ITOPJ) GO TO 115 

KSTART = ITOPI 

IF II .EO. KSTART) GO 70 125 
IF (KSTART .GT. KEND) GO TO 125 
LKI = LITOPI-1 
LKJ = L ITCP J+ITOP I— ITOPJ— 1 
GO TO 12C 
115 KSTART = ITOPJ 

IF (I .EO. KSTART) GO TO 125 
IF (KSTART .GT. KEND) GO TO 125 
LKI = L ITCP1+ITCP J-I TOPI— 1 
LKJ = LITOPJ-1 
120 DO 122 K=KSTART,KEND 
LKI = LKI +1 
LKJ = LKJ+1 

122 S = S— V (K)*V(LKI ) *V( LK J ) 

125 V(LIJ) = S 

IF (I .LE. ND) V(LIJ) = S/VII) 

128 LITOPI = LITPPI+LVIl) 

134 IF (J .CT. ND) GO TO 140 
V( J) = V( LJJ) 

IF (ITOPJ .EQ. J) GO TO 129 
LKJ = LITOPJ-1 
DO 138 K= ITOPJ » JM 1 
LKJ = LKJ +1 

138 V(J) = V(J) - V(K )*V(LKJ)**2 

139 NERR0R=4 
IF (ABSIV(J) ) .LT • EPS ) GO TO 999 

V(LJJ) = 1.0 

140 CONTINUE 

C GROUP A OPERATE ON GPOUP P. 

C I COLUMNS ARE IN GROUP A, J COLUMNS IN GROUP L. 

IF (IGA .EC. NG) GO TO 195 

IGAP1 = IGA+1 

JEGP = JEGA 

DO 19? IGP = IGAP1 ,NG 

NCGR = LV (K V02+IGB) 

JSGP = JEGF + 1 
JEGB = JSGB+NCGB— 1 
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LEGB = LSGB-1 
DO 1 S J= JSGBt JEGP 
151 LFGB = LFGB+LV( J ) 

CALI YIN (NUTP,V,LSGB,LEGB) 

LJJ = LSGP-1 

DO 190 J= JSGBt JEGB 

JM1 = J-l 

ITOPJ = J-LV(J)+I 

LITOPJ = LJJ + 1 

LJJ = LTTOPJ+LV ( J )— 1 

IF ( ITOPJ .GT. JEGA) GO TO 190 

I ST APT = ITOPJ 

LIJ = LITOPJ-1 

IF (ITOPJ .GF. JSGA) GO TO 155 
ISTART = JSGA 
LIJ = LIT OP J ♦ JSGA— ITOPJ— 1 
155 LITOPT = LSGA 

IF (ISTART .FO. JSGA) GO TO 160 
ISM1 = ISTART— 1 
DO 157 1= JSGA» ISM! 

157 LITOPI = LITOPI+LV(I) 

160 DO 178 I=ISTART, JEGA 
KFND -1~1 
IF (I .GT. ND) KEND=ND 
LIJ = LIJ+1 
S = V(LIJ) 

1 Ml = 1-1 

I TOP 1 = I -LV ( I ) + 1 

IF (ITPPI .LT. ITOPJ) GO TO 165 

KSTART = ITCPI 

IF (I .FC. KSTART) GO TO 175 
IF (KSTART .GT. KEND) GO TO 175 
LKI = L ITPPI— 1 
LKJ = LITOPJ+nOPl-ITOPJ-1 
GO TO 170 
165 KSTART = ITOPJ 

IF (I .FO. KSTART) GO TO 175 
IF (KS1AP1 .GT. KFND) GO TO 175 
LKI = L ITPPI +IT0PJ-1T OP 1—1 
LKJ = LITCPJ-1 
170 DP 17 2 K=KSTART » KEND 
LKI = LKI +1 
LKJ = LKJ +1 

172 5 = S— V (K )*V ( LKI T*V ( LK J ) 

175 V(LIJ) = S 

IF (I .LE. ND) V ( LIJ) =S/V( I) 

178 LITOPI = LITOPI+LV( I ) 

190 CONTINUF 

19? CALL YOUT (NUTQ ,V,L SGP , LFGB ) 

195 CALL YOUT (NUT3 »V* LSGA »LFGA ) 

RFWIND NUTA 

CALL Y0UT1 iNUTA»LV • 1 »KV) 

C CONVERT (U) FROM HAND (NUT3) TP SPARSE (NUT1 ) NOTATION. 

C DISASSEMBLE TO GET (022) ON NUTR. 
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REWIND NUT3 
REWIND NtfTI 
LVGS = KV -NGRCUP 
LVR = LVGS 

DC 202 1GR0UP=1 » NGROUP 
LVR = LVR +1 

202 LV(LVP) = LV IKV02+IGR0UP ) 

LS = LVGS-NRA 
LVE = LS 
DO 204 1= 1 ,NRA 
LVE = LVF+1 
204 LV(LVF) = LV ( I ) 

KVMAX = KV/4 

JF (KVMAX. GT. IS) KVMA X=LS 

MHFAD(l) = NRA 

MHF AD ( 2 ) - NPA 

MHEAD13 ) = NGRCUP 

MHEAD (4 ) = NNZZ 

MHEADC5) = 0 

MHFAD C 6 ) = 0 

MHF AD ( 7 ) = 5HUPPER 

CALL YOUT t (NUT1 ,MHE AD, 1* 10 ) 

LVT = C 
LVR = LVGS 
LVF = LS 
LVFP = LS 
IZ - 0 

DC 250 IGR0UP=1, NGRCUP 
LVR = LVR +1 
LZ = 0 

NROWS = LV(LVR) 

NELG = C 

DO 206 1*1, NROWS 
LVE * LVF + 1 

206 NELG = NFLG+LV( LVF) 

CALL YIN (NUT3,V,1,NELG) 

DO 20P 1=1, NROWS 
IZ = IZ + 1 
LVEP = LVFP+1 
JS = IZ— L V( LVFP ) +1 
DO 208 JZ=JS,1Z 
LZ = LZ+1 

208 LV C LZ ) = 1C0CC0*JZ+IZ 
MPHEAD(l) = LZ 
MPHFAD ( 2 I = LVU ) 

MPHFAD ( 3) •- LV(LZ) 

CALL YOUTI (NUT1 ,MPHEAD,1 ,10) 

CALL YOUT I (NUT1 ,LV, 1 *LZ) 

CALL YOUT (NUT 1 ,V, 1 , LZ ) 1 

250 CONTINUE 

CALL YPAFT (NVTl ,V,L V,K V,NUT ? i 

CALL YD IS A INUT1 ,ND+ 1 ,ND+1 ,NUTR »NR ,NR,V,LV,KV,NUT2) f 

t CALCULATE REDUCTION TRANSFORMATION MATRIX. 

IF ( IFT .ED. t) RETURN 
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C ( U ) IS ON NUT3 • BANDED, NO MATRIX OR PARTITION HEADERS. 

: TRANSFER Ull TO NUTT (TOP NON-ZERO IN COLUMN DOWN TO DIAGONAL) AND 

C U12 TO NIH1 (INDIVIDUAL FULL COLUMN). 

REWIND NUT4 

CALL YIM (NUT4,LV,I,KV) 

do 302 :-i,ic 

302 MPHEAD(l) = 0 
REWIND NUT3 
REWIND NUTT 
REWIND NUT1 
JUE = 0 
NGU11 = 0 
DO 335 IG = 1 ,NG 
NCG = LV( KV02+IG ) 

JUS s JUE +1 
JUE = JUS+NCG-1 
NFLG = 0 

DO 306 JU=JUS, JUE 

306 NELG = NELG + LV(JU) 

CALL YIN (NUT 3 »V, 1 »NFLG ) 

IF 1ND .GF. JUE) CO TO 310 

IF (ND .GF. JUS) GO T 0 320 

IF (ND .LT. JUS) GO TO 330 

NERROR = 5 
GO TO 999 

310 MPHEAD ( 1 ) = JUS 
MPHEAD ( 2 ) = JUF 
MPHE AD ( 3 ) = NELG 
CALL Y0UT1 (NUTT,MPHEAD,1 ,10) 

CALL YOUT (NUTT ,V,1 , NELG) 

NGU11 = NGUll+1 
GO TO 335 
320 NFL ' 0 

DO 322 JU= JUS »ND 

322 NEL = NFL + LV(JU) 

MPHEAD ( 1 ) = JUS 
MPHEAD (2) = ND 
MPHEAD ( 3 ) = NEL 

CALL YPUTI (NUTT »MPHFAD,1 ,1.0) 

CALL YOUT (NUTT ,V,1 , NEL) 

NGU11 = NGU11+1 

IF (JUS . FU . JUE) GO TO 335 

LFJll = NFL 

JUSX = ND+1 

323 DO 327 JU=JUSX,JUF 
ITOP = JU-t V ( JU ) ♦ I 
LSJU = LEJU+1 

LFJU = LS JU+LV ( JU )-l 
DO 3?A IV- 1 , ITOP 
LIV = NFLG+IV 

324 V(LIV) - 0.0 

IF (ITOP .GT. ND) GO TO 327 
LUV * LSJU-1 
DO 325 1 V=ITOP ,ND 
LIV = NELG+IV 
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LUV = LUV+1 
325 V(LIV) = V(LUV) 

327 CALL YOU! (NUT1 ,V,N ELG+1 ,NELG+ND ) ] 

GO TO 335 
330 LF JIJ = 0 
JUSX = JUS 
GO TO 323 
335 CONTINUE 

C GROUP (U12) INDIVIDUAL FULL COLUMNS FROM NUT1 ONTO NUT2. 

C USE VCI THPU (KV-N1/2) TO AGREE WITH YBSL3A. 

RFWIND NUT1 
REWIND NUT? 

K VI - (KV-N)/2 
LF = 0 
NCG = 0 
NGP = 0 
00 343 J=1,NR 
LS - LF+1 
LE = LS+ND-1 

CALL YIN (NUTI,V,LS,LE) 

NCG = NCG+1 

IF (J .EO. NR) GO TO 342 
IF ((LF+ND) .LF. KV1 ) GO TO 343 

342 MPHFAD ( 1 ) = NCG 

CALL YOUTI (NUT2,MPHEAD,1,10) 

CALL YOUT (NUT2,V,1,LE) 2 

LE --- 0 

NCG = 0 

NGP = NGB + 1 

343 CONTINUF 
C 

C BACK SOLUTION FOR (T) FROM (Ull )*(T)=(UI2) . 

C fUU) GROUPS ARE OBTAINED IN BACKWARDS ORDER. 

C V(1 THRU (KV-NJ/2 ) CONTAINS Y=U12, 2=1 COLUMNS OF A GROUP. 

C V((KV-NI/?+l THRU KV-N) CONTAINS COLUMNS OF U (FROM TOP NON— ZERO 
C DOWN TO DIAGONAL) OF A GROUP. 

C LV( I ) » 1=1 »N IS NUMBER OF ELEMENTS IN COLUMN I. 

LSU = (KV-N)/2 + 1 
REWIND NUT2 
REWIND NUT 1 
C 

DO 389 IGR=1,NGP 

CALL YINI (NUT 2 »IH » 1 . 1C ) 

NCIGP = I H ( 1 ) 

NFLIGB = ND*NCIGB 

CALL YIN (NUT? f Vtl tNELlGB ) 

DO 357 1GUX=1 ,NGU11 

BACKSPACE NUTT 

BACKSPACE NUTT 

CALL YINI (NUTT ,IH, 1,1C) 

JSU = Ih(l) 

JEU = 1H ( 2 ) 

NFL IGU = IH(3 ) 

CALL YIN (NUTT »V,LSU,LSU+NF LIGU—1 ) 

BACKSPACE NUTT 
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BACKSPACE NUTT 

DO 35? JP =1 »NC IGF 

L2SM1 = (JR-1)*ND 

LITJU = LSU+NELIGU 

DO 356 JUX=JSU,JEU 

JU = JSU+JEU-JUX 

LJJU = LITJU-1 

LITJU = L JJU-LVC JU) + 1 

IT JU = JU-LVt JUJ+1 

IF (1TJU .EQ. JU) GO TO 356 

LJJUM1 = LJJU-1 

12 = LZSMT+JU 

L2Y = L2SM1+ITJU— 3 

DO 354 LU=LITJU, L JJUM 1 

LZY = LZY*" 

354 V ( L ZV I = V(LZY) - V(LU)*V(LZJ 

356 CONTINUE 

357 CONTINUE 

DO 359 1GU=1 ,NGU 1 1 

CALL YINI (NUTT ,IDUM ,1 , 1 ) 

359 CALL YIN (NUTT, DUM,l,l) 

DC 372 1=1,10 
372 IH ( I ) = 0 

IH(1 ) = NCIGE 

CALL Y0UT1 (NUT1 »IH» 1,10) 

DO 375 1=1 ,NEL IGB 
375 V(I) = -V(I) 

CALL YOUT ( NUT 1 *V , 1 » NE Ll G6 ) 1 

389 CONTINUE 

CONVERT (T) FROM FULL COLUMN (NUT1) TO SPARSE (NUT3) NOTATION. 

REWIND NUTT 

REWIND NUT 3 

IH ( 1 ) = ND 

IH ( 2 ) = NR 

IH( 3 ) = NGfl 

1 H (4 ) a ND*NR 

IH( 5 ) = C 

IH ( 6 ) = C 

I H ( ? ) = 5HWHCLE 

CALL YOUTI (NUT3 t lH*I ,10) 

JZ = 0 

DO 395 IGP=1,NGB 

CALL Y1N1 (NUT 1 ,IH, 1,10) 

NC = IH ( 1 ) 

NN2PB = ND*NC 

CALL YIN (NUT1 ,V,1 ,NN2PB ) 

LB = G 

DO 392 J =1 ,NC 
JZ = JZ+1 
DO 392 I Z = 1 , ND 
LB * LB +1 

392 LV (LB) = 1COOGO*I2+JZ 
IH ( 1 ) = NNZPB 
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IH ( 2 ) = LVH ) 

IHI3) = LVfNNZPP) 

CALL YPtm CNUT3 *IH* 1,10) 

CALL VPUT1 (NU^ 3 ,LV, 1 ,NNZPB ) 

395 CALL YOUT INUT3, V,1,NNZPP) 3= 

CALL YNOZER INUT3 ,V,LV,KV,NUT1 ) 

CALL YLPPD <NUT3,V,LV,KV,NUT1,NUT2) 

CALL YZEPO INUTT # NPA,NR ) 

CALL YASSEM <NUT3, 1 ,1 ,NUTT, V,LV,KV,NUT1 ,NUT2,NUT4) 

CALL YUNITY (NUT3 ,NR, V,LV,KV) 3= 

CALL YASSEM <NUT3,ND+1 ,1 ,NUTT, V.LV,KV,NUT1 ,NUT2,NUT4) T= 

RETURN 

999 CALL ZZE Mb I6HYSRED2 ,NERROR ) 

END 
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SUBROUTINE YSTCD (NUT A , A ,NRA ,NCA, KR A, KCA, V *LV ,KV, NUT1 ) 

DIMENSION V(1)»LV(1),A(KRA*1) * MHE AD (10) ,MPHEaDI10) 

DATA NIT, NOT/5, 6/ 

C 

C CONVERT A MATRIX FROM SPARSE NOTATION TO DENSE NOTATION. 

C CALLS FORMA SUBROUTINES YIN ,YINI ,ZZPOMB. 

C DEVELOPED BY R A PhILTPPUS. JANUARY 1<»60. 

C LAST REVISION BY WA BENFIELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS 

C NUTA = INPUT LOGICAL NUMBER OF UTILITY TAPE ON WHICH SPARSE MATRIX 
C A IS STORED. 

C A = OUTPUT DFNSE MATRIX. SIZE (NRA *NCA ) • 

C NRA - OUTPUT NUMBFR OF ROWS IN A. 

C NCA = OUTPUT NUMBER OF COLUMNS IN A. 

C KRA = INPUT ROW DIMENSION OF A IN CALLING PROGRAM. 

C KCA = INPUT COLUMN DIMENSION OF A IN CALLING PROGRAM. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV = INPUT DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = INPUT LOGICAL NUMBER OF UTILITY TAPE . 

C 

C NERROR EXPLANATION 

C 1 = DIMENSION SIZE EXCEEDED (KRA, KCA). 

C 2 = DIMENSION SIZE EXCEEDED (KVI. 

C 3 = ROW OR COLUMN LOCATION OUTSIDE MATRIX A. 

h 4 = ROW OR COLUMN LOCATION OUTSIDE MATRIX A. 

REWIND NUT A 

CALL YINI (NUTA ,MH£ AO, 1 , 10 ) 

NRA = MHEAD(l) 

NCA = MHE AD (2) 

NERR0R=1 

IF (NRA.GT.KRA .OR. NCA.GT.KCA) GC TO 999 
C 

DO 10 1=1, NRA 
DO 10 J=I,NCA 
10 A( I , J)=0. 

C 

NPART = MHE AD( 3 ) 

1 SHAPE = MhE AD ( 7 ) 

DO 40 K-I , NPART 

CALL YINI (NUT A ,MPHEAD ,1,10) 

NNZP = MPHEAD(l) 

NERROR =2 

IF (NNZP.GT.KV) GO 10 999 
IF (NNZP.G7.0) GO TO 20 
CALL YINI (NUT A »LV ,1,1) 

CALL YIN (NUT A ,V , 1 , 1 ) 

GO TO 4C 
C 

20 CALL YIM (NUT A ,LV, 1 ,NNZP ) 

CALL YIN (NUTA,V,1 ,NNZP) 

DO 30 L-l ,NNZP 
I=LV(L )/10CC00 
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J=LVCL)-10C000*I 

NERR0R=3 

IF II.GT.KRA .HR. J.GT.KCA) GC TO 999 
IF C ISHAPE *EO. 5HWFOLE ) GO TO 3C 

NF.RROR=4 

IF CJ.GT.KRA .OR. I.GT.KCA) GO TO 999 
A ( Jy I ) =V( L ) 

30 A Cl » J)=V( L) 

40 CONTINUE 
RETURN 

999 CALL ZZBOMB ( 5HYS10D ,NERRCR) 

END 
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SUBROUTINE YSYMLH (NUTAZ »V*LV,KV*NUT1 ,NUT2 ) 

DIMENSION Vfll ,LV<1),MHEA0(10) 

DATA NIT, NOT/5, 6/ 

C 

C SYMMETRIZE SPARSE MATRIX AZ BY PLACING VALUES FROM ABOVE THE DIAGONAL 
C BELOW THE DIAGONAL. 

C CALLS FORMA SUBROUTINES YIN ,Y1NI ,YLORD ,YNOZER ,YOLT ,YOUTI , 

C YPART . 

C DEVELOPED BY R A PHILIPPUS. JUNE 1969. 

C LAST REVISION BY R A PHILIPPUS. JUNE 1973. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTAZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX AZ IS STORED. 

C V = VECTOR WORK SPACE. 

C LV = VECTOR WORK SPACE. 

C KV - DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

C NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

C NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

C 

CALL YPART (NUTAZ, V, LV,KV,NUT1 ) 

C 

REWIND NUTAZ 
REWIND NUT1 

CALL Y INI (NUT AZ*MH EAD,I ,10 J 
NNZA = MHEAD(A) 

IF (NNZA.FQ.C) RETURN 
ISHAP = MHEAD(T) 

IF ( ISHAP. E0.AHD1AG) RETURN 
MHEAD(7> = 5HWH0LE 
C 

NPAR » m - MHEAD (3 ) 

DO 20 1=1 »NPARTA 

CALL YINI (NUT AZ,MHEAD ,1C ,101 

NNZP = MHFAD(IO) 

CALL YINI (NUT AZ, LV, 1 , NNZP) 

CALL YIN (NUTAZtVs I, NNZP) 

NNZQ=NNZP 

C 

DO 10 J=1 *NNZP 
I A=LV ( J )/lCOOOC 
J A=l V ( J )- ICOGGO* I A 
IF (Ia.EC.Ja) GO TO 10 
IF UA.C-T.JA) GO TO 5 
NNZ0=NNZ0+1 
NNZ A=NN2 A+l 
L V ( NNZ Q )= I CCOCC*JA+IA 
V (NN20 ) =V ( J ) 

GO TO 10 
5 V(J)=C. 

10 CONTINUE 
C 

MHEAD(IO) = NNZO 

CALL Y0U7I (NUT) , MHEAD, 1C, 10) 

CALL YOUTI (NUT 1 ,LV, 1 , NNZO) 

20 CALL TOUT (NUT1 ,V, I ,NNZQ ) 
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RFWIND NUTAZ 
RFWIND MUT1 
MHFAD(4> = NNZA 
MHFADC5) = 0 
MHFAD (6 ) = 0 
MHEAD(IO) = 0 

CALL YOUTI (NUTAZ, MHEAD, 1,10 
C 

DO 25 1=4,10 
25 MHEAP(I) = 0 

DO 30 1=1 ,NPARTA 

CALL YIN1 (NUT1 ,MHE AD , 1 , 1 ) 

NN2.P = MHFAD (1) 

C ALL YINI (NUT1 ,LV, 1 ,NNZ P) 

CALL YIN (NU11 ,V,1,NNZP) 

MHEAD ( 2 ) = LV(1) 

MHFAD ( 3 ) = LV(NNZP) 

CALL YOUTI (NUTAZ, MHEAD, 1,10) 

CALL YOUTI (NUTAZ, LV»1,NNZP) 

3C CALL YOUT (NUTAZ, V, 1 ,NNZP ) 

C 

CALL YNOZER (NUTAZ, V, LV,KV,NUT1 ) 

CALL YLORD ( NUTAZ, V, LV,KV»NUT1,NUT?) 

RFTURN 

END 
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C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 


c 


c 


SUBROUTINE YSYMUH (NUTAZ ,V,LV,KV,NUT1 ,NUT2 ) 

DIMENSION V(1),LV(1), MHEAD ( 10 ) 

DATA NIT ,NCT/5 *6/ 

SYMMFTR IZE SPARSE MATRIX AZ BY PLACING VALUES PROM BELOW THE DIAGONAL 
ABOVE THE DIAGONAL. 

CALLS FORMA SUBROUTINES YIN ,YIN1 * YLORD ,YNOZER,YOUT ,YOUTI , 

YPART . 

DEVELOPED BY R A PHILIPPUS. JUNE 1969. 

LAST REVISION EY R A PHILIPPUS. JUNE 1973. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTAZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX AZ IS STORED. 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

CALL YPART (NUTAZ, V * LV,KV»NUT1 ) 

REWIND NUTAZ 
PEWIND NUT1 

CALL YINI (NUTAZ, MHEAD, 1*1G) 

NNZA = MHEAD(4) 

IF (NNZA.FO.C) RETURN 
ISHAP = MHEAD(7) 

IF (ISHAP. E0.4HDIAG) RETURN 
MHEAD (7) =• 5.HWHOLE 

NPARTA = MHEAD ( 3 ) 

DO 20 1=1, NPARTA 

CALL YINI ( NUT AZ, MHEAD, 1C, 10) 

NNZP = MHEAD CIO) 

CALL YINI (NUTAZ, LV,1, NNZP) 

CALL YIN (NUTAZ»V , 1 ,NNZP ) 

MNZQ=NNZP 

DO 10 J=1,NNZP 
IA=LV( J)/100000 
JA=L V( J )- 100000* I A 
IF (IA.FO.JA) GO TO 10 
IF (IA.LT.JA) GO TO 5 
NNZQ=NNZQ+1 
NNZ A=NNZ.A+1 
LV(NNZO )=ICOCOO*JA+IA 
V (NNZQ ) =V ( J ) 

GO TO 10 
5 V(J)=0. 

10 CONTINUE 

MHEAD(IO) = NNZO 

CALL Y0UT1 (NUT 1 ,MH EaD ♦ 1 C , 1 C ) 

CALL YOUTI (NUTI ,LV,1,NNZQ) 

20 CALL YOUT (NUTI ,V,1,NNZQ) 
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$ 

REWIND NUTAZ 
RFWIND NUT1 
MHEA0(4) = NNZA 
MHEAD (5 ) = 0 
MHEAP ( 6 ) = 0 
MHEAP(IO) = 0 

CALL YOUTI ( NUT A 2., MHEAD ,1,10) 

C 

DO 25 1=4,10 
25 MHEAD(I) = 0 

DO 30 1 = 1 ,NPART A 

CALL YINI (NUT 1 ,MHE AD, 1 » 1 ) 

NNZP * MHEAD(l) 

CALL YINI (NUT1 ,LV,1,NNZP) 

CALL YIN (NUT1 ,V,1 ,NN2P) 

MHEAD (2 ) = LV(1) 

MHEAD( 3 ) = LV(MNZP) 

CALL YOUTI (NUTAZ, MHEAD, 1,10 
CALL YOUTI (NUTAZ, LV,1*NNZP) 

30 CALL YOUT (NUTA Z,V , 1 ,NNZP ) 

C 

CALL YNCZER ( NUTAZ, V , LV,KV, NUT 1 ) 

CALL YLORO (NUTAZ,V,LV,KV,NUT1,NUT2) 

RETURN 

END 
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SUBROUTINE YTRANS (NUT A,NUTAT , V, L V, KV ,NU7 I ,NUT2 ) 

DIMENSION V(l) ,LV(1) » MHEAD (10) ,MPHEAD(10) 

DATA NIT »NCT/5 ,6/ 

C 

C TRANSPOSE SPARSE MATRIX A INTO SPARSE MATRIX AT. 

C CALLS FORMA SUBROUTINES YIN ,YINI ,YLORD ,VCUI ,YOUTI ,YPART , 

C ZZBOMB. 

C DEVELOPED BY R A PHILIPPUS. JANUARY 1960. 

C LAST REVISION BY WA PENFIELD. MARCH 1976. 

C 

C SUBROUTINE ARGUMENTS tALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORFD. 

NUTAT = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX AT IS STORED. 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION S ! Z c OF V,LV IN CALLING PROGRAM. 

NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

NUT2 =• LOGICAL NUMBER OF UTILITY TAPE. 

NERROP EXPLANATION 
1 = DIMENSION SIZE EXCEEDED (KV). 

REWIND NUTA 
REWIND NUTAT 

CALL YINI (NUT A ,MHEAD» 1,10) 

N = MHFAD(l) 

MHEAD ( 1 ) = MHEAD ( 2) 

MHEAD (2 ) = N 
I ASHA P = MHEADI7) 

1 SHAP=I ASHAP 

IF ( IASHAP.EC.5HUPPFR ) I SHAP=5HLCWER 
IF (IASHAP.EC.5HL0WER) 1SHAP=5HUPPER 
MHEAD ( 5 ) = 0 
MHEAD ( 6 ) = 0 
MHFAD ( 7 ) = ISHAP 
CALL YOUTI (NUTAT,MHEAD,1,10) 

NNZ A a MHEAD (A ) 

NPART = MHE AD ( 3) 

IF (NNZA.GT.O) GO TO 3 
DO 7 1=1, ir 
7 MPHFAD(l) = 0 

CALL YOU T I (NUTAT,MPHEADf I ,10) 

CALL YOUTI (NUT AT, MP HEAD ,1,2) 

CALL YOUTI (NUTAT, MPHEAD, 1,2) 

RETURN 
C 

3 DO 10 1=1, NPART 

CALL YINI (NUTA,MPHEAD,1 ,10) 

NNZP = MPHEAD(l) 

NERRCR=1 

IF (NNZP.GT.KV) GO 70 999 
CALL YINI (NUTA, LV, 1, NNZP) 

CALL YIN (NUTA »V» 1 ,NNZP ) 

C 

DO l J=1 ,NNZP 



YTRANS — 2/ 2 


IA=LV( J )/ 100000 

5 LVC J)=1C0000*(LV( J)-1C0000*IA ) +1 A 
C 

CALL YCUTI (NUT AT*MPHEAD»1»10) 

CALL YOUTI (NUTAT.LV, 1,NNZP ) 

CALL YCUT (NUTAT,Vt l.NNZP) 

10 CONTINUE 
C 

CALL YLORD <NUTAT,V t LV,KV,NU1 1,NUT2) 
RETURN 
C 

999 CALL ZZB0M6 (6HYTRANS * NERROR ) 

END 
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SUBROUTINE YUNITY (NUTA, NRA, V,LV,KV) 

DIMENSION V(1),LV(1),MHEAD(10) 

GENERATE SPARSE UNITY MATRIX A. CONES ON THE DIAGONAL) « 

CALLS FORMA SUBROUTINES YOUT ,YOUTJ . 

DEVELOPED BY R A PHIL1PPUS. JANUARY 1<>70. 

LAST REVISION BY R A PHIL1PPUS. AUGUST 1973. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH UNITY MATRIX A IS 
STORED. 

NRA = SIZE OF UNITY MATRIX A ( SQUARE ) • 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

REWIND NUTA 

NPARTA= (NRA-D/CKV/A) +1 

LVC I ) = NRA 
LVC2 ) = NRA 
LVC3 ) = NPARTA 
LVC A) = NRA 
LVC 5 ) = 5HORDER 
LVC6 ) = KV 
LVC7) = AHDIAG 
LVC8 ) = 0 
LVC9) = 0 
LVC 10) = C 

CALL YOUTT (NUTA .LV, 1 , 10 ) 

DO 5 1= A, 1C 
5 MHEADC I) = 0 
LAE=KV/A 
J=0 

DO 10 1=1, NRA 
J=J+1 

LVC J ) =I 00000*1-*- 1 
V C J ) - 1 . 

IF CJ.LT.LAE .AND. l.LT.NRA) GO 10 10 
MHEADC 1) = J 
MHEADC?) = LVC 1 ) 

MHEADC3) = LVC J) 

CALL YOUTI (NUT A ,MHE AD , ) , 10 ) 

CALL YOUTI (NUT A ,LV» 1 ,J) 

CALL YOUT (NUT A ,V, 1 , J ) 

J = 0 

10 CONTINUE 
C 

RETURN 

END 
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SUBROUTINE YWRITE (NU TA, AN AME , V, LV , KV ) 

DIMENSION V(1),LV(1),W(1G) ,MHE AD( 1C ) 

DATA NIT , NOT/5 ,6/ 

C 

C WRITE SPARSE MATRIX A ON PAPER IN SAME FORMAT AS DENSE FORMA 
C SUBROUTINE WRITE. RECUIRES 132 COLUMN (MINIMUM) PRINTER. 

C UP TO 1C DATA FIELDS PER LINE. PRINT ONLY NON— ZERO FIELD ROWS. 

C CALLS FORMA SUBROUTINES PAGEHD, YIN ,YINI , ZZBOMB . 

C DEVELOPED BY R A PHILIPPUS. SEPTEMBER 196E. 

C LAST REVISION BY WA BENFIELD FOR NASA. MAY 1976. 

C 

C SUBROUTINE ARGUMENTS (ALL INPUT) 

C NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX A IS STORED 
C ANAME = MATRIX IDENTIFICATION. < A6 FORMAT) 

C V = VECTOR WORK SPACE. 

C LV - VECTOR WORK SPACE. 

C KV = DIMENSION SIZE OF VyLV IN CALLING PROGRAM. 

C 

C NEPROR EXPLANATION 

C 1 ~ DIMENSION SIZE EXCEEDED (KV). 

C 

2001 FORMAT (//2?H OUTPUT SPARSE MATRIX A6 ,2XIH ( 15 ,2H XI5,6H ) ( 15 

* 2 1H N0N-2F.R0 ELEMFNTS) (I5.12H PARTITIONS ) 1XA6// 

* 10X,10(7XylH(yl2,1H) )/) 

2002 FORMAT (//22H OUTPUT SPARSE MATRIX A6 ,2X1H ( 15 ,2H XI5,6H ) ( 15 

* 21H NON-ZERO ELEMENTS) (I5,12H PARTITIONS ) 1XA6, 10H CONTINUED// 

* 10X,10(7X,1H( ,I2ylH) )/) 

2003 FORMAT ( 1 X, 2 15 ,2 X ,1P1 0E1 1 .4) 

2004 FORMAT (15H0END OF YWRITE.) 

3001 FORMAT (45H0END OF YWRITE. NRA OR NCA HAS BEEN EXCEEDED 15, 

* 7H TIMES.) 

PULL UP A NEW PAGE FOR MATRIX AND PRINT MATRIX NAME. 

REWIND NUTA 

CALL YINI (MUTAyMHEADyl ylO) 

NRA = MHFAD(l) 

NCA = MHFAP ( 2 ) 

NPART = MMFAD(3) 

NNZA = MHEAD (4) 

MCKORD = MHE AD ( 5 ) 

KVCHK = MHEAD ( 6 ) 

ISHAPF MHFAD(7 ) 

CALL PACEHD 

WRITF (NOT, 2001) *NAM E , NRA , NCA , NNZA , NPART , 1SH APE ,( 1, 1=1 , 10 ) 

IF ( NNZA • EC .0 ) GO TO 40 
NL INF = 0 
I FL AC-0 
I JK=0 
C 

DO 3B M^lyNPARI 

CALL YINI (NUTA, MHEAD, 1,10) 

NNZP = MHFADH ) 

LFFLP = MHEAD ( 2 J 
LLELP = MHF AD ( 3 ) 

IF ( NNZP.GT . 0 ) GO TO 2 
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CALL Y1NI (NUTA »MHE AD » 10 ,10 ) 

CALL YIN (NUTA, V,KV,KV) 

MHFADflC) = 0 
GO TO 38 

2 NERROR=l 

IF (NNZP.GT.KV) GO TO 999 
Cl 4.L YINI CNUTA.LV, I ,NNZP) 

CALL YIN (NUTA,V,1,NNZP) 

c 

DO 35 1=1 ,NNZP 
I A=LVC T )/ 100000 
JA-LV(I)— 1 00000*1 A 

IF (IA.GT.NRA .OR. JA.GT.NCA) IJK=IJK+1 
IF (I.EC.l .AND. M.EQ.l) GO TO 20 
K= JA-JS+1 

IF (IA.NE.IS .OR. K.LE.O .OR. K.GT.10) GO TO 5 
W(K) = Vm 

IF (I.LT.NWZP .OR. M.LT.NPART) GO TO 35 
IFLAG=1 
5 N J=1 0 

IF C (JS+9).GT.NCA) NJ=NCA-JS+1 
IF (JA.GT.NCA) N J=10 
NLINF=NLINF+1 
IF INLINE. LE.4M GO TO 10 
CALL PAGE HD 

WRITE (N0T,2CC2) ANAME .NRA.NCA.NNZA .NPART , 1 SHAPE » (J»J=1»10) 
NLINF = 1 

10 WRITE ( NOT . 2003 ) IS , J S , < W ( J) , J=1 ,N J ) 

IF (NNZP.FQ.C) GO TO 38 
IF (IFLAG.EO.l) GO TO 35 
C 

C SKIP A SPACE BETWEEN EACH ROW IF THERE ARE MOPE THAN 1C COLUMNS 
C AND SOMETHING HAS BEEN WRITTEN. 

IF (NCA.LF.10 .OR. IL.EG.IA) GO TO 25 
NLINE = NLINF + 1 
WRITE (NOT , 2C03 ) 

20 IL=IA 
25 I S=I A 

JS=( JA-I )/I0*10+l 
DC 30 L = 1 ,10 
30 W ( L ) =0. 

K=JA-JS+1 

W(K)=V(I) 

IF (I.LT.NNZP .OR. M.LT.NPART) GO 10 35 
IFLAG=1 
GO TO 5 
35 CONTINtJF 
38 CONTINUE 
C 

IF (IJK.EO.O) GO TO 40 
WRITE (NOT, 3001) 1JK 
RETURN 

40 WRITF ( NOT, 20C4 ) 

RETURN 



999 CALL ZZBOMB <6HYWR1TE »NERROR ) 
END 
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SUBROUTINE YWTAPE (NUTA, ANAME * V,LV,EV,NTAPE) 

DIMENSION V( I)»LV(1)»MCHECM2),MHEAP(10) 

COMMON /L START/ IRUNO t DATE ,NPAGE ,UN AME { 3) » TITLE! ( 12) »TITLfc2 (12) 
DATA NIT tNOT/5 *6/ 

OATA BUF » EOT ,NGNE *SPART / 0. , 3HE0T, 1 , 5HSPART / 

WRITE SPARSE MATRIX AON TAPT (NTAPE). 

INITIALIZE TAPE WITH SUBROUTINE INTAPE . 

REWIND TAPE BEFORE F1PST USE OF THIS SUBROUTINE. 

NOTE. ..THIS ROUTINE IS DESIGNED SPECIFICALLY FOR WRITING ON A DISK 

(EG CDC-64C0 DISK). USING THIS ROUTINE 70 WRITE ON A PHYSICAL 
TAPE DIRECTLY (IE WITHOUT USING THE DISK AS AN INTERMEDIARY) 
WILL PROBABLY GIVF POOR RESULTS ( DUF TO THE TOLERANCE 
CHARACTERISTICS OF A TAPE DRIVE) AND SHOULD BE AVOIDED IF AT 
ALL POSSIBLE. 

...THE CDC-B4CC DISK IS AU'f 0MA7 ICAl LY ENDFILED AFTER FACH WRlit. 
CALLS FORMA SUBROUTINES YIN ,Y1NI ,Z2B0MB. 

DEVELOPED BY R A PHILIPPUS. NOVEMBER ! c f>8. 

LAST REVISION BY Wa 6ENF1ELD. MARCH 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA = LOGICAL NUMBER OF UTILITY TAPF ON WHICH MATRIX A IS STORED. 
ANAME = MATRIX IDENTIFICATION. ( A6 FORMAT) 

V = VECTOR WORK SPACE. 

LV = VECTOR WORK SPACE. 

KV " DIMENSION S } 2E Of 1 V,LV IN CALLING PPOGRAM. 

NTAPE * LOGICAL NUMBER. OF TAPE ON WHICH MATRIX A IS TO RE WRITTEN. 

NERROP EXPLANATION 
1 - DIMENSION SIZE EXCEEDFD (KV). 

INTERNAL VARIABLES THAT ARE PUT ON TAPE (TRANSFERRED THRU COMMON). 
RUNNO IS PUN NUMBER OF PROBLEM. ( A6 FORMAT). 

DATE IS DATE. ( A6 FORMAT). FOR EXAMPLE ISFEfcS 

SEARCH TAPF FOR END OF WRITTEN DATA. 

10 READ (NTAPE) TAPE ID,LN ,1 EOTCK 
IF ( IEOTC K. EC.3HE OT)GO TO 20 
READ (NTAPE) 

GO TO 10 

END OF WRITTEN DATA HAS BEEN FOUND. 

20 BACKSPACE NTAPE 
REWIND NUTA 

CALL YINI (NUTA,M, HEAD, 1,10) 

NR A = MHFAD(l) 

NCA -- MHF AD ( 2 ) 

NPART = MHF AD ( 3 ) 

NN2 A a MHEAD ( 4 ) 

MCHECK C! ) = MHEAD (5 ) 

MCHFCK ( 2 ) s MHEAD (6) 

MSHAPE = MH F A D ( 7 ) 

IF (NPART. GT.O) GO TO 25 

WRITE (NTAPE ) TAPFID, LN.BUF, IRUNO, ANAMF ,NR A, NCA, DATE. ,SPART,feUF r 

* NONE, NONE tMSHAPEt (BUF, '•=1 ,6) 
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WRITF (NT APE ) BUF ,BUF 

LN=LN*1 

GO TO 4C 

25 DO 35 J-l tNPAPT 

CALL Y1NI <NUTA»NHEAD* It 10 ) 

NN2P = MBEAD(l) 

NERRC»R=1 

J* (NN2P.GT.KV) GO TO °99 

WRITE (NT APE 1 TAP c IO, LNtBUF,IRUNOtANAMEtNRA,NCA,DATF,SPAPT ,NNZP»J» 
* NP ART » ( MC HECK ( I ) ,1=1 ,2) ,MSHAPE, (FUFtI-lt4> 

IF (NNZP.GT.O) GO TO 30 
CALL YJNI (NUT A »NHE AD*10,10) 

CALL YINI (NUTA, V, 1 , 1) 

NHFAD(IC) = C 
WRITE (NT APE) EU F ,BUF 
GO TO 35 

30 CALL VIM (NL'T A ,LV , 1 ,NNZP1 
CALL Yiw (NL’TA,V,1,NNZP) 

WR ITE (NT APE) (LV (I ) , V C 1 ) , 1 = 1 ,NN2P) 

35 LN=LN+1 

40 WRITE (NTAPf) TA PE1D, LN ,ECT, (BUF*I=1 , 16 I 
BACKSPACE NTAPE 

return 

999 CmLL ZZECMB (6HYV/TAPE ,NEF.RCSR ) 

END 
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SUBROUTINE YZERLH (NUTAZ ,V,LV,KV,NUT1 ,N'UT2 ) 

DIMENSION VII ) ,L V (1 ) , MHE AD 1 1C » 

DATA \’IT,NCT/5,6/ 

ZERO THF LOWFR HALF CF SPARSE MATRIX AZ. 

CALLS FORMA SUBROUTINES YIN ,YINI , YLORD ,YCUT ,YCUTI ,YPART , 

ZZBOMB. 

DEVELOPED BY R A PHILIPPUS. JUNE I«69. 

LAST REVISION BY WA BENFIELD FOR NASA. NAY 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTAZ = LOGICAL NUMPEP OF UTILITY TAPE ON WHICH MATRIX AZ IS STORED. 

V = VECTOR WORK SPACE. 

LV = VECTOR WC'PK SPACE. 

KV = DIMENSION SIZE OF V,LV IN CALLING PROGRAM. 

NUT1 = LOGICAL NUMBER OF UTILITY TAPE. 

NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

NERRCF EXPLANATION 

1 = MATRIX NOT SOUAPF. 

2 = DIMENSION SIZE EXCEEDED (KV). 

GET (A) HEADER INFORMATION. 

REWIND NUTAZ 
REWIND NUT 1 

CALL YIN1 (NUT AZ,MHE AD *1 ,1 0 ) 

NPA = MHEAD(l) 

NCA = MHEAD ( 2 ) 

NERROR=I 

IF (NRA.NE.NCA) GO TO 999 
NPA»7A = MHEAD (3 ) 

NN2. A = MHEAD (M 
MCKORD = MHF AD ( 5 ) 

IASHAP = MHE AD ( 7 ) 

IF ( NN Z A .EO. 0) RETURN 
IF (IASHAP. EC.^HDIAG) RFTURN 
IF ( IASHAP. EC-.5HUPPER ) RETURN 
MSHAPE = 5HUPPFR 

IF (IASHAP. FO.BHLCWER) MSHAPE=4HDI AG 

NNZZ=0 

NPARTZ=0 

LOOP C.' ; (A) PARTITIONS. 

DO 20 1 = 1 ,NPARTA 

CALL YINI (NUTAZ, MHEAD, 1 ,10) 

NNZP = MHEAD ( I ) 

NERR0R=2 

IF (NNZP.G7.KV) GO TO 999 
CALL YTNI (NUT AZ,L V , 1 ,NNZP ) 

CALL YIN (NUT AZ, V, 1 ,NNZP ) 

DC 10 J=I ,NNZP 
IA-LV( J)/1C000G 
J A=LV ( J ) — I PC COO* I A 
IE (JA.LT.IA) V(J)=0. 
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10 CONTINUE 
NNZQ=C 

DO 15 J=1,NNZP 
IF (V(J).FC.C.) GC TO 15 
NNZQ=NNZO+T 
VI NNZO ) =V ( J ) 

LV(NNZQ)=LV(J> 

15 CONTINUE 

IF (NNZC.FO.C) GO TO 20 

NNZZ=NNZZ+NN20 

MHEADtl) = NNZO 

CALL YCUTI (NUT1 ,MHEAD* 1 » 1 ) 

CALL YOUTI (NUT1 »LV» 1 ,NNZG ) 

CALL YOUT (NUT 1 »V» 1 »NNZQ ) 
NPARTZ=NPARTZ+1 
20 CONTINUE 

TRANSFER OATA FROM NUT1 TO NUTAZ. 

IF (NNZZ.FC.NNZA) GO TO 40 

RFWINO NUTAZ 

REWIND NUT1 

MHF AD ( 1 ) = NRA 

MHEADC?) = NCA 

MHE AD ( 3 ) = NPARTZ 

MHEADfM = NNZ2 

MHFAD ( 5 ) = MCKCRD 

MHE AD ( 7 ) = MSHAPE 

CALL YOUTI (NUTAZ, MHEAD,I, 10) 

DO 3C 1=1, NPARTZ 

CALL YINI ( NUT 1,MHEAD,I,1) 

NNZO = MHEAD(l) 

CALL YINI (NUT I »LV,1,NNZC) 

CALL YIN (NUTT ,V,1 , NNZO) 

MHFAD (2 ) = LV(1J 
MHFAD(? ) = LV(NNZQ) 

MHF AD (4 ) = 0 
MHFAD ( 5 ) = 0 
MHE AD ( 7 ) = 0 

CALL YOUTI (NUTAZ, MHFAD, 1 ,10) 

CALL YOUTI (NUTAZ, LV,1, NNZO) 

30 CALL YOUT (NUTAZ, V, 1 , NNZO ) 

40 CALL YLORD (NUTAZ, V, LV,KV,NUT1,NUT2 ) 
RETURN 

999 CALL ZZB0M6 (6HYZFRLH ,NERROR ) 

END 
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SUBROUTINE YZERO (NUTA ,NR A,NCA ) 

DIMENSION MHEADliO) 

GENERATE A NULL SPARSE MATRIX A. 

CALLS FORMA SUBROUTINES YOUT ,YOUTI . 

DEVELOPEO BY R A PH1LIPPUS. OCTOBER 1969. 

LAST REVISION BY JOHN ADMIRE *NASA* EEB 1974. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTA = LOGICAL NUMBER OF UTILITY TAPE ON WHICH NULL MATRIX A IS 
STORED. 

NR A = MUMBFP OF ROWS IN A. 

NCA = NUMBER OF COLUMNS IN A. 

REWIND NUTA 
MHEAD(I) = NR A 
MHEAD ( 2 ) = NCA 
MHEAD (3 ) = 0 

MHEAD (4) = 0 

MHEAD ( 5 ) = 0 
MHEAD (6) = 0 
MHF AD 1 7 ) = BHWHOLE 
MHEAD(B) = 0 
MHEAD ( Q ) = 0 
MHF AD ( 1 0 ) = 0 

CALL YOUTI { NUT A »MHE AD » 1 » 1C ) 

MHEAD(l) = C 
MHE AD ( 2 ) " 0 
MHEAD ( 7 ) = 0 

CALL voUTT (NUT A »MHEAD» 1* 10 ) 

RETURN 

END 
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SUBROUTINE YZERUH (NUTAZ ,V ,LV ,KV, NUTI ,NUT2 ) 

DIMENSION V(I)»LV(1)»MHEAD(10) 

DATA NIT, NOT/5, 6/ 

ZERO THE UPPER HALF OF SPARSE MATRIX AZ. 

CALLS FORMA SUBROUTINES YIN ,Y1NI , YLORD ,YOUT ,YGUT1 ,YPART , 

ZZBOMB. 

DEVFLOPED BY R A PHILIPPUS. JUNE 1969. 

LAST REVISION EY WA BENFIELD FOR NASA. MAY 1976. 

SUBROUTINE ARGUMENTS (ALL INPUT) 

NUTAZ = LOGICAL NUMBER OF UTILITY TAPE ON WHICH MATRIX AZ IS STORED. 

V = WORK VECTOR . 

LV = WORK VECTOR. 

KV = DIMENSION SIZE OF V,LV IN CALL INC PROGRAM. 

NUTI - LOGICAL NUMBER OF UTILITY TAPE. 

NUT2 = LOGICAL NUMBER OF UTILITY TAPE. 

NERPQR EXPLANATIONS 

1 = MATRIX NOT SC'UARF. 

2 = DIMENSION SIZF EXCEEDED (KV). 

GET (A) HEADER INFORMATION. 

REWIND NUTAZ 
REWIND NUTI 

CALL Y1N1 (NUTAZ, MHEAD, 1,10) 

NR A = MHEAD(l) 

NCA = MHE AD ( 2 ) 

NERRGR=1 

IE (NPA.NF.NCA) GO TO 999 
NPARTA = MHE AD ( 3 ) 

NNZ A = MHFAD(A) 

MCKORD = MHE AD ( 5 ) 

IASHAP = MHE AD (7) 

IF (NNZ A .EO. 0) FETURN 
I F (IASHAP. EQ.4HD1AG) RFTUPN 
IF ( IASHAP. EC. EHLCWER ) RETURN 
MSHAPE = 5HL0WER 

IF ( IASHAP. EQ.BHUPPER ) MSHAPE = 4HDIAG 

NNZZ=C 

NPARTZ=C 


LOOP ON (A) PARTITIONS. 

DO ?C 1 = 1, NPARTA 

CALL YINI (NUTAZ, MHEAD,1 ,10) 

NNZP = MHEAD(l) 

NERRCR=2 

IF (NNZP. GT . KV ) GO TO 999 
CALL YINI (NUTAZ, LV,1, NNZP) 

CALL YIN (NUT AZ,V , X ,NN2P ) 

DO 10 J=I,NNZP 
I A=LV ( J )/ 100000 
JA=LV( J)-I0C0CP*IA 
IF (JA.GT.IA) V(J)=0. 
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10 CONTINUE 
NNZ0=0 

DC 15 J=1,NNZP 

I F (V(J).EO.O.) GO TO 15 

NN20-NN7C' - * - ! 

V I NNZO ) =V ( J ) 

L V ( NNZO )=LV ( J ) 

15 CONTINUE 
C 

IE (NNZO. EC. C) GC TC 2 C 

NNZZ=NNZZ+NNZO 

MHEAD) 1 ) = NNZO 

CALL YCUTI (NUT1 ,MHF AD * 1 , 1 ) 

CALL YCUTI (NUT1 ,LV,1,NNZ0) 

CALL YCUT (NUT1 ,V, 1 ,NNZQ ) 
NPARTZ=NPART2+1 
20 CONTINUE 
C 

C TRANSFER DATA FROM NUT1 TO NUTAZ. 

IF (NNZ2.FC.NN2A) GO TO AC 

REWIND NUTAZ 

REWIND NUT1 

MHE AD ( 1 ) = NRA 

MHE AD ( 2 ) = NCA 

MHEAD(?) = NPARTZ 

MHE AD ( A ) = NNZZ 

MHE AD (5 ) = MCKORD 

MHEAD(6) = C 

MHE AD ( 7 ) = M SHAPE 

CALL YCUTI (NUTAZ,MHEAD,1, 1C) 

C 

DO 3G 1=1, NPARTZ 

CALL YTNT (NUT1 ,MHE AD, 1, 1 ) 

NNZO = MHEADd) 

CALL YINI (NUT) »LV* 1 » N N 2 0 ) 

CALL YIN (NL'Tl ,V,1 ,NNZC ) 

MHEAD ( 2 1 = LV(1) 

MHEAD ( ? ) = LV(NNZQ) 

MHFAD(A) = C 

MHEADd) = 0 
MHEADJ7) = C 

CALL YCUTI (NUTAZ, MHEAD, ) ,)C ) 

CALL YCUTI ( NUT AZ,LV,1, NNZO) 

30 CALL YCUT (NUTAZ, V, 1 , NNZO ) 

C 

AO CALL YLORD (NUTAZ, V , LV,KV , NUT 1 ,N’JT2 ) 
RETURN 
C 

999 CALL ZZBOMF ( 6HYZ ERUH , NERRuP ) 

END 



